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1.  Introduction 

This  report  describes  the  mathematical  analysis  on 
which  the  program  modifications  to  Fotonap  are  based.  Also 
included  in  the  report  is  some  program  documentation  and  a 
description  of  the  final  test  runs  used  in  the  program  check 
out.  The  program  modifications  may  be  subdivided  into  7 separate 
groups,  6 of  which  are  mathematical  in  nature  and  therefore 
described  here.  To  deal  with  the  added  capabilities,  the  set 
of  control  cards  used  to  run  Fotonap  has  been  augmented.  A new 
Fotonap  user's  guide  has  been  issued  as  a separate  document. 
Program  changes  based  on  the  analysis  given  in  this  report  have 
been  implemented  and  checked  out  on  both  the  CDC  6400  and  the 
Uni vac  1108  versions  of  Fotonap. 

The  most  fundamental  change  to  Fotonap  is  the  inclusion 
of  a Kalman  filter  and  a fixed  lag  smoother.  The  smoother 
formulation  is  considerably  more  complicated  than  that  of  the 
filter.  This  applies  both  to  the  mathematical  analysis  and  to 
the  program.  The  core  storage  requirements  are  also  much  higher. 
The  analysis  is  based  on  Gelb  (1977),  though  checking  back  to 
the  original  source  (Medltch,  1969)  one  of  the  required  equations 
was  found  to  be  in  error.  The  formulation  implemented  in  Fotonap, 
however,  differs  slightly  from  that  given  by  Meditch  (1969). 

The  Fotonap  formulation,  though  mathematically  equivalent, 
requires  slightly  less  core  storage.  The  filter  and  smoother 
analysis  is  given  in  Section  2.  The  description  of  the  satellite 
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drag  coefficient  and  the  GPS  user  clock  as  random  variables  is 
based  on  some  equations  given  by  Dr.  Ballew  (1977)  of  DMA  Aero- 
space Center,  St.  Louis.  This  is  given  in  Sections  3 and  4. 

Formulae  for  calculating  the  state  transition  matrices 
(required  by  the  filter/smoother)  from  those  computed  by  the 
regular  Fotonap  integrator  are  given  in  Section  5. 

A brief  description  of  the  Global  Positioning  System 
(GPS)  is  given  in  Section  6.  A derivation  of  all  the  equations 
required  for  handling  GPS  measurements  is  also  contained  within 
this  section.  A considerable  amount  of  effort  was  spent  on 
developing  an  efficient  scheme  for  selecting  GPS  satellites  when 
running  the  program  in  a simulation  mode.  It  had  originally  been 
thought  possible  to  use  the  Morduch  (1976)  method,  but  that  was 

developed  for  somewhat  different  requirements.  The  rather  lengthy 
analysis  is  given  in  Appendix  E. 

Section  7 gives  the  mathematical  analysis  for  drag 
segmentation,  which  is  a scheme  whereby  the  satellite  drag 
coefficient  changes  discretely  at  fixed  (by  the  program  user) 
intervals  of  time.  A novel  (optional)  feature  of  the  implemented 
scheme  is  that  the  coefficients  may  be  mutually  constrained,  the 
strength  of  the  constraints  being  determined  by  the  program  user. 
(Very  strong  constraints  effectively  eliminate  the  distinction 
between  the  segments.  This  feature  was  utilized  in  the  program 
check-out. ) 

A new  atmospheric  model,  the  Lockheed- Jacchia  Atmosphere 
has  been  added  to  Fotonap.  The  program  is  a modification  of  one 
supplied  by  Mr.  George  Stentz  (1978)  of  DMA  Aerospace  Center.  A 


J 

description  of  the  formulae  used  is  given  in  Section  8.  Section 
9 of  the  report  gives  a general  description  of  all  program  changes. 

The  changes' to  the  Fotonap  user’s  guide  are  indicated  in  Section 
10.  Section  11,  the  last  section  of  the  main  text,  describes 
some  of  the  test  runs  made  in  checking  out  the  new  version  of 
Fotonap. 

A set  of  7 appendices  are  also  included  in  the  report. 

The  first  5 (A  through  E)  give  detailed  derivations  of  some  of 
the  formulae  used  in  the  main  text.  Appendices  F and  G updates 
some  of  the  previous  program  documentation  (Hartwell,  1975) 
relating  to  Fotonap. 

Fotonap  is  sometimes  spelled  Fhotonap  in  this  report. 

It  is  the  same  program. 
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Mathematical  Analysis  For  Fixed  Lag  Smoother 


Definition  of  Terns 

E Operator  denoting  'Expected  value  of' 

A-1  Inverse  of  matrix  A 
T 

A Transpose  of  matrix  A 

A‘T  = (A"1)1  - (A*1)  - 1 

Cov  y Covariance  of  y 

T 

Cov  y = E(y-Ey) (y-Ey) 

x(k,n)  Estimate  of  parameter  vector  at  time-point 
k after  processing  measurements  at  time- 
points  1 through  n. 

xl(k)  = x(k,k) 

x2(k)  = x(k+l,k) 

xS(k)  = x(l,k)  if  k s L 

• x(k-L,k)  if  k > L 

L Lag  constant 

b(k,n)  Estimate  of  parameter  vector  at  time-point 
k after  processing  measurements  at  time- 
points  n through  N and  n £ k 

xT(k)  True  parameter  vector  at  time-point  k 

x(k,n)  = x(k,n)-xT(k) 

b(k,n)  - b(k,n)-xT(k) 

P(k,n)  Estimated  covariance  of  x(k,n) 

Pl(k)  - P(k,k) 

P2(k)  - P(k+l,k) 


PS(k) 


P(l,k)  if  k < L 


P(k-L.k) 


if  k > L 
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B(k,n)  Estimated  covariance  of  b(k,n) 

M(j,k)  State  transition  matrix  relating  the 

parameter  vectors  (system  states)  at  time- 
points  j and  k. 

M(k,j)  = M(j.k)-1 

W(k)  State  noise.  See  equations  (2.4.1)  and  (2.4.2) 

Q(k)  = E W(k)W(k)T 

A(k)  = P(k,k)M(k+l,k)'rP(k+l,k)-1 

C(j  ,k)  = A(j  )A(j+l)  . . . . A(k)  [Defined  only  for  k > j] 

Cl(k)  = C(h,k),  where  h=l  (mod  L) 

and  k-L<h<k 

C2(k)  = C(k+1-L,k)  [Defined  only  for  k > L] 

dP(k)  Change  in  covariance  of  parameter  estimate 
at  time-point  k [See  Equations  (2.2.6)  and 
(2. 2. 7)1 

dx(k)  Change  in  parameter  estimate  at  time-point 
k [See  equation  (2.2.13)  and  2.2.14)] 

AP(k)  = P(k-L.k)  - P(k-L,k-L) 

= PS(k)  - Pl(k-L) 

Ax(k)  = x(k-L,k)  - x(k-L,k-L) 

= xS(k)  - xl(k-L) 

G(k)  Kalman  gain  matrix  [See  equation  (2.2.8)] 

H(k)  Matrix  of  partial  derivatives  of  measurements 

with  respect  to  parameters 

Z(k)  Vector  of  measurements 

r Measurement  noise 

T 

R = Err 

F(k)  = P(k,N)  - P(k,N+l) 
f(k)  = x(k,N+l)  - x(k,n) 
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2. 2 Filter  and  Smoother  Equations 

The  following  formulae  will  be  derived. 

P (k+1 , k)  = M(k+1 ,k)P(k,k)M(k+l,k)T+Q(k) , (2.2.1) 

A(k)  = P(k,k)M(k+l,k)TP(k+l,k)“ 1 (2.2.2) 

C(j ,k)  = C(k,k-l)A(k)  (2.2.3) 

C(j  ,k-l)  = AO-D-'CO-l.k-l)  (2.2.4) 

C(j.j-l)  = I.  the  identity  matrix  (2.2.5) 

P (k+1, k+1)  “ P(k+l,k)  - dP(k+l)  (2.2.6) 

dP(k+l)  - G (k+1) H (k+1) P (k+1 , k)  (2.2.7) 

G(k+1)  = P (k+1 , k) H (k+1) T [R(k+1)  + 

+ H(k+l)P(k+l,k)H(k+l)Tf  * (2.2.8) 

P(l,k+1)  = P(l,k)  - C(l,k)dP(k+l)C(l,k)T  (2.2.9) 


P(k+1-L,k+1)  - P(k+1-L,k-L)  + 

+ A(k-L)_1AP(k)A(k-L)-T 

- C(k+l-L,k)dP(k+l)C(k+l-L,k)T  (2.2.10) 

AP(k)  - P(k-L,k)  - P(k-L.k-L)  (2.2.11) 


1 


i 
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x(k+l,k)  = M(k+l,k)x(k,k)  (2.2.12) 

x(k+l,k+l)  = x(k+l,k)  + dx(k+l)  (2.2.13) 

dx  (k+1 ) = G(k+l)[Z(k+l)-H(k+l)x(k+l,k)]  (2.2.14) 

x(l,k+l)  = x(l ,k)  + C(l,k)  dx (k+1)  (2.2.15) 

x(k+l-L,k+l)  = x(k+l-L,k-L)  4-  A(k-L)- 1 Ax(k) 

+ C(k+1-L,k)  dx(k+l)  (2.2.16) 

Ax(k)  = x(k-L.k)  - x(k-L,k-L)  (2.2.17) 


The  significance  of  the  most  important  of  the  above 
equations  may  be  described  as  follows : 

Equations  (2.2.1)  and  (2.2.12).  Covariance  and  parameter 
vector  propagation  from  time-point  k+1  in  the  absence  of  any 
measurements  passed  time-point  k. 

Equations  (2.2.6)  and  (2.2.13).  Covariance  and  parameter 
vector  update  at  time-point  k+1,  after  the  measurements  at  time- 
point  k+1  have  been  processed. 

Equations  (2.2.9)  and  (2.2.15).  [Fixed  point  smoother] 
Estimated  covariance  and  vector  at  the  initial  time-point  after 
the  measurements  at  the  first  k+1  time-points  have  been  processed. 
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Equations  (2.2.10)  and  (2.2.16).  [Fixed  lag  smoother] 
Estimated  covariance  and  vector  L time  intervals  prior  to  latest 
measurement  point.  As  can  be  seen  from  equation  (2.2.11)  and 
(2.2.17)  the  smoother  estimate  for  the  previous  time-point  is 
always  needed  in  the  calculations.  The  fixed  point  smoother 
is  used  to  get  such  an  estimate  for  the  initial  time-point. 


| The  following  definitions  are  introduced 


Pl(k)  = P(k,k)  (2.3.1) 

P2(k)  = P (k+1 , k)  (2.3.2) 

PS(k)  - P(l,k)  if  k < L (2.3.3) 

PS(k)  = P(k-L,k)  if  k > L (2.3.4) 

Cl(k)  =»  C(h,k),  where  (2.3.5) 

h = 1 (mod  L)  (2.3.6) 

and 

k-L  < h < k (2.3.7) 

C2(k)  - C(k+1-L,k)  (2.3.8) 

xl(k)  = x(k,k)  (2.3.9) 

x2(k)  - x(k+l,k)  (2.3.10) 

xS(k)  = x(l,k)  if  k < L (2.3.11) 

xS(k)  = x(k-L,k)  if  k > L (2.3.12) 


Equations  (2.2.1)  through  (2.2.17)  may  now  be  rewritten 
in  the  new  notation.  Since  there  is  no  risk  of  confusion,  the 
indices  for  Q,  M,  G,  H,  R,  dP,  AP,  dx  and  Ax  will  be  dropped.  We 
find  that 

P2(k)  = M Pl(k)MT  + Q (2.3.13) 

A(k)  = Pl(k)MTP2(k)-*  (2.3.14) 

Cl(k)  - A(k)  for  k = 1 (mod  L)  (2.3.15) 

Cl(k)  = Cl(k-l)A(k)  for  k jf  1 (mod  L)  (2.3.16) 
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C2(k)  = A(k-L)- 1 C2(k-l)A(k)  for  k f 0 (mod  L)  (2.3.17) 


C2(k)  = Cl(k)  for  k = 0 (mod  L) 


(2.3.18) 


[Note  that  repeated  use  of  equation  (2.3.17)  for  the 
computation  of  C2  will  result  in  numerical  inaccuracy. 
C2  is,  therefore,  reset  every  L cycles  using  equation 
(2.3.18)] 


Pl(k+1)  = P2(k)  - dP 


(2.3.19) 


dP  = GH  P2(k) 


(2.3.20) 


G = PHT[R  + H P2(k)  HT]_1 


- C2(k)  dP  C2(k) 


(2.3.21) 


PS(k+l)  - PS(k)  - Cl(k)  dP  Cl(k) 1 f or  k < L (2.3.22) 


PS(k+l)  - P2(k-L)  + A(k-L)_1AP  A(k-L)"T 


f or  k > L (2.3.23) 


AP  = PS(k)  - Pl(k-L) 


(2.3.24) 


x2(k)  =M  xl(k) 


(2.3.25) 


xl(k+l)  = x2(k)  + dx 


(2.3.26) 
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dx  = G lZ  - H x2(k)] 

(2.3.27) 

xS(k+l)  = xS(k)  + Cl(k)dx, 

for  k < L 

(2.3.28) 

xS (k+1)  - x2(k-L)  + A(k-L) “ 

*Ax  + C2(k)dx 

• 

for  k > L 

(2.3.29) 

Ax  = xS(k)  - xl(k-L) 

(2.3.30) 

It  can  be  seen  from  the  above  that  variables  Pi,  P2, 


A,  xl  and  x2  require  L storage  blocks  each,  but  for  the  remaining 
quantities  only  the  last  computed  value  need  be  maintained  in 
computer  memory. 


Derivation  of  equations 


The  state- transition  matrix  M(k+l,k)  relates  the 
true  parameter  vector  xT(k)  at  time-point  k to  the  corresponding 
vector  at  time-point  k+1  through  the  equation 

xT(k+l)  - M(k+l,k)xT(k)  + W(k),  (2.4.1) 

where  the  term  W(k)  arises  due  to  our  lack  of  knowledge  of  the 
system.  W(k)  is  thus  unknown  to  us.  However,  we  shall  assume 
that 

EW(k)  = 0,  EW(k)W(k)T  = Q(k)  (2.4.2) 

and  W(k)  and  W(m)  are  assumed  to  be  uncorrelated  if  k f m.  For 
the  forward  prediction  formula  we  obtain  a similar  formula 

x(k+l,k)  = M(k+l,k)x(k,k)  (2.4.3) 

[Same  as  equation  2.2.12.] 
whence  the  error  x must  satisfy, 

x(k+l,k)  = M(k+l,k)x(k,k)  - W(k)  (2.4.4) 

The  covariance  propagation  equation  (2.2.1)  follows  from  the 
above  and  equation  (2.4.2). 


The  backward  prediction  formula  corresponding  to  equation 


(2.4.3)  is  given  by 

b(k,k+l)  = M (k, k+1) b (k+1, k+1)  (2.4.5) 

with  the  error  b being  given  by 

b(k,k+l)  « M(k,k+1)  [b(k+l,k+l)  + W(k)]  (2.4.6) 

From  the  above  and  equation  (2.4.2)  it  follows  that  the  covariance 
is  given  by 

B(k,k+1)  = M(k,k+l)[B(k+l,k+l)  + Q(k)] M(k,k+1)T  (2.4.7) 

whence 

B(k+1 , k+1)  = M (k+1 , k) B (k , k+l)M (k+1 , k) T - Q(k)  (2.4.8) 

Adding  equations  (2.2.1)  and  (2.4.8)  we  obtain 

P (k+1 , k)+B (k+1 , k+1)  « M (k+1 , k) [ P (k , k)+B (k , k+1) ] M (k+1 , k)T 

(2.4.9) 

Since  the  forward  parameter  estimate,  x(k,k),  is  based 
on  measurements  1 through  k,  and  the  backward  estimate,  b(k,k+l), 
is  based  on  measurements  k+1  through  N,  it  follows  that  the  two 
estimates  are  statistically  independent.  We  may  thus  use  formulae 
(B.l)  and  (B.2)  in  Appendix  B to  get  the  smoothed  solution  at 
time-point  k: 


A 
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x(k,N)  = P(k,N) [P(k,k)-1x(k,k)+B(k,k+l)-1b(k,k+l)]  , (2.4.10) 

P(k,N)  = [P(k,k)“ *+  BCk.k+l)"1]"1  (2.4.11) 


The  above  argument  also  applies  to  estimates  x(k,k-l) 
and  b(k,k),  which  may  be  combined  to  yield  formulae  similar  to 
(2.4.10)  and  (2.4.11).  At  time-point  k+1  we  thus  obtain 


x(k+l,N)  = P(k+1,N) [P(k+l,k)"1x(k+l,k)+B(k+l,k+l)“1 

b (k+1, k+1)],  (2.4.12) 

P(k+1,N)  = [P(k4-l,k)"1+B(k+l,k+l)_1J"1  (2.4.13) 

From  the  above  equation  we  deduce  that 

P(k+1,N)  = B (k+1, k+1) [P (k+1, k)  + B(k+l,k+l)]  *P(k+l,k) 

= P (k+1 , k)  - P (k+1 , k) [ P (k+1 , k)  + B(k+l,k+l)l -IP(fc+l,k) 


From  the  above  and  equation  (2.4.9)  we  obtain 


P(k+1,N)  = P (k+1 , k) 

-P(k+l,k)M(k+l,k)'T[P(k,k)  + B(k,k+1)] 
M(k+l,k)“1P(k+l,k) 

Since  by  equation  (2.2.2) 

P (k+1 , k)M(k+l , k) “T  = A(k)_1P(k,k). 


(2.4.14) 


(2.4.15) 
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i£Su‘! 


it  follows  that 

P(k+1,N)  = P(k+l,k) 

-A(k)_,P(k.k)[P(k,k)  4-  B(k>k+l)]‘1p(k,k)A(k)"T  (2.4.16) 

From  equation  (2.4.11)  we  find  that 

P(k,N)  - P(k,k)[  P(k,k)  + B(k,k+1)  ] " ^(k.k+l) 

- P(k,k)  - P(k,k)[P(k,k)  4-B(k,k4-l)  ] “ ’p(k.k) 

(2.4.17) 


From  the  above  and  equation  (2.4.16)  we  deduce  that 

P(k+1,N)  = P(k+l,k)  4-  A(k)“1lP(k,N)-P(k,k))A(k)"T 

(2.4.18) 

Let 

F (k)  - P(k,N)  - P (k,N4-l)  (2.4.19) 

It  then  follows  from  equation  (2.4.18)  that 

F (k4-l)  - A(k) " 1 F(k)A(k)_T  (2.4.20) 

and  hence  that 

F (k)  - A(k)  F(k4-1)  A(k)T  (2.4.21) 

From  the  above  we  deduce  that 

F (k)  - A(k)A(k4-l).  . A(N)F(N+1)A(N)T.  . .A(k)T  (2.4.22) 

For  j<k,  Define  C(j,k)  by 

C(j.k)  = A(j)A(j+l).  . .A(k)  (2.4.23) 
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It  can  easily  be  seen  that  the  above  definition  is  consistent 
with  equations  (2.2.3)  through  (2.2.5).  From  equations  (2. A. 22) 
and  (2.4.23)  we  obtain 

F (k)  = C(k,N)F(N+l)  C(k,N)T  (2.4.24) 

It  follows  from  equations  (2.4.19)  and  (2.4.24)  that 

P(k,N+l)  - P(k,N)-C(k,N)[P(N+l,N)-P(N+l,N+l)] 

C(k,N)T  (2.4.25) 

Replacing  k by  1 and  N by  k in  the  above  equation  we  obtain 

P(l,k+1)  - P(l,k)  - C (1 ,k)dP (k+l)C(l ,k)T  (2.4.26) 

where 

dP(k+l)  = P(k+l,k)  - P(k+l,k+l)  (2.4.27) 

The  derivation  of  equations  (2.2.6)  through  (2.2.8)  is  given  in 
Appendix  C.  Since  equations  (2.4.27)  and  (2.2.6)  are  equivalent 
it  follows  that  equations  (2.4.26)  and  (2.2.9)  are  identical. 

In  order  to  derive  equation  (2.2.10)  we  first  substitute 
k-L  for  k and  k for  N in  equation  (2.4.18)  thus  obtaining 

P(k-l-L,k)=P(k+l-L,k-L)+A(k-L)- 1 AP(k)A(k-L)“T  (2.4.28) 

where  AP(k)  is  given  by  equation  (2.2.11) 


A 
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In  equation  (2.4.25)  we  substitute  k+l-L  for  k and  k 

for  N: 

P(k+l-L,k+l)=P  (k+1- L, k) -C (k+l-L, k)dP  (k+1)  C (k+l-L,k)T  (2.4.29) 

where  SP(k+l)  is  defined  by  equation  (2.4.27).  Equations 
(2.4.28)  and  (2.4.29)  may  be  seen  to  be  equivalent  to  equation 
(2.2.10).  This  completes  the  derivation  of  equations  (2.2.1) 
through  (2.2.12).  Equations  (2.2.13)  and  (2.2.14)  are  derived  in 
Appendix  C. 

We  shall  now  proceed  to  derive  equations  (2.2.15)  through 
(2.2.17).  It  follows  from  equations  (2.4.12)  and  (2.4.13)  that 

x(k+l,N)  - b (k+1, k+1) 

+ P (k+1 , N) P (k+1 , k) “ 1 [ x (k+1 , k) -b (k+1 , k+1)]  (2.4.30) 

From  the  above  and  equations  (2.4.18),  (2.4.3)  and  (2.4.5)  we 
deduce  that 

x(k+l,N)  - b (k+1, k+1)  + [x (k+1 , k) -b (k+1 , k+l)J 
+A(k)-1[P(k,N)-P(k,k)]  A(k)-Tp(k+l,k)"1  • 

M(k+1 , k) [ x (k , k) -b (k , k+1)  ] (2.4.31) 

whence  by  equation  (2.2.2), 

x(k+l,N)  = x(k+l,k) 

+A(k)~  1 [P(k,N)-P(k,k)j  P (k,k) “ 1 [x(k,k)-b (k,k+l)J, 
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x(k+l,N)  « x(k+l,k)+A(k)“ '{ -x(k,k)  + 

P(k,N) {P(k,k)~ 1 [x(k,k)-b(k,k+l)]  + P(k,N)‘  ^(k.k+l) }}, 

whence  by  equations  (2.4.11)  and  (2.4.10) 


x(k+l,N)  - x(k+l,k)+A(k)"1[x(k,N)-x(k,k)]  (2.4.32) 

Let  f (k)  = x(k,N+l)  - x(k,N)  (2.4.33) 

Hence  f(k+l)  = A(k)“*f(k)  (2.4.34) 

and  f (k)  = A(k)  A(k+1) . . . A(N)f (N+l)  (2.4.35) 

Using  equation  (2.4.23)  we  deduce  that 

f (k)  = C(k,N)f (N+l)  (2.4.36) 

Hence , 

x(k,N+l)  = x(k,N)  + C(k,N)dx(N+l) , (2.4.37) 

where 

dx(N+l)  » x (N+l, N+l)  - x(N+l,N)  (2.4.38) 


If  in  the  above  two  equations  we  substitute  1 for  k and 
k for  N we  obtain  equation  (2.2.15)  and  the  equivalent  of  equation 
(2.2.13).  In  order  to  derive  equation  (2.2.16)  we  first  substitute 
k-L  for  k and  k for  N in  equation  (2.4.32),  thus  obtaining 


x(k+l-L,k)  = x(k+l-L,k~L) 

+A(k+L)“ 1 Ax(k) , (2.4.39) 

where  Ax(k)  is  given  by  equation  (2.2.17).  In  equation  (2.4.37) 
we  substitute  k+l-L  for  k and  k for  N: 

x(k+l-L,k+l)  = x (k+l-L, k)  + C (k+l-L, k) dx (k+1) , (2.4.40) 

where  dx(k+l)  is  defined  consistently  with  equation  (2.2.13). 

It  is  easily  seen  that  equations  (2.4.39)  and  (2.4.40)  may  be 
combined  to  yield  equation  (2.2.16).  This  completes  the  derivation 
of  the  filter  equations. 
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3.  Modelling  of  the  Drag  Coefficient 


The  drag  coefficient  is  modelled  as  a constant  plus 
a small  perturbation  xp,  which  satisfies  the  differential  equation 

xp(t)  = -xp(t)/Tp  + vp(t),  (3.1) 

where  a dot  denotes  differentiation  with  respect  to  t and  vp 
is  a random  variable  with 

E vp(t)  = 0,  and  (3.2) 

E vp(tl)vp(t2)  = 6(t2  - tl)  2Qp/Tp,  (3.3) 

where  6(t)  is  the  Dirac  delta  function. 

Making  the  substitution 

s = t/Tp,  (3.4) 


and  defining 


xp(t)  = x'  ( s)  , (3.5) 

where  a prime  denotes  differentiation  with  respect  to  s,  we 
deduce  that 


x"(s)  + x'(s)  = Tp  vp(t)  (3.6) 

If  further  we  define  v(s)  = Tp  vp(t)  then  it  follows  from 
equations  (3.2),  (3. 3)  and  equation  (D.38)  in  Appendix  D that 
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E v(s)  = 0,  and 


(3.8) 


Ev(sa)v(sb)  = 6 [(sa-sb)Tp]Tp  2Q/Tp 

= 2Q  6 (sa-sb)  (3.9) 

Comparing  equations  (3.5)  through  (3.9)  with  equations 
(D. 1)  through  (D. 3)  we  deduce  from  equations  (D.7)  and  (D. 10) 
that 


xp(t)  = gp  xp(t0)  + hp ( t ) , (3.10) 

where 

gp  = exp[-(t-t0)/Tp]  , (3.11) 

Ehp(t)  = 0,  (3.12) 

and 

E hp (t) 2 = Qp [l  - gp2]  (3.13) 


1 


Equations  (3.10)  through  (3.13)  may  be  seen  to  correspond 
to  equations  (2.4.1)  and  (2.4.2).  In  order  to  interpret  Qp  and 
Tp  we  note  that  if  t is  much  greater  than  to  then  gp  is  negligibly 
small.  We  hence  find  from  equations  (3.10)  and  (3.13)  that 

E xp ( t ) 2 = Qp  (3. 14) 

From  equation  (D.33)  we  deduce  with  the  aid  of  equation  (3.4)  and 
(3.5)  that 
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5.*p(t)_xp<t+at)  . exp(_4t/Tp)  (3. 15) 

E xp(t)xp(t) 

This  concludes  Che  description  of  the  modelling  of  the  drag 
coefficient. 

3. 1 Simulation  of  Drag  Coefficients. 

Equations  (3.10)  through  (3.13)  are  used  in  the  simulation 
of  drag  coefficients,  fip(t)  is  chosen  as  a normally  distributed 
random  variable  satisfying  equations  (3.12)  and  (3.13). 


i 
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4.  Modelling  of  Clock  Bias 

The  clock  bias  b(t)  is  modeled  as  the  integral  of  the 
clock  bias  rate  b(t),  which  itself  satisfies  the  differential 
equation 

b(t)  = - b(t)/Tb  + vb(t),  (4.1) 

where  a dot  denotes  differentiation  with  respect  to  t and  vb  is 
a random  variable  with 


E vb(t)  = 0,  and  (4.2) 

E vb(t  ,)vb(t  2)  = 6(t2  - t ,) 2 Qb/Tb2  (4.3) 

where  6(t)  is  the  Dirac  delta  function. 

Making  the  substitution 

s - t/Tb  (4.4) 

defining 

b(t)  = x(s)  (4.5) 

and  denoting  differentiation  with  respect  to  s by  a prime,  it 
follows  that 

x"(s)  + x’(s)  - v (s ) , (4.6) 

where 

v(s)  - vb ( t)  Tb2  (4.7) 


-23- 


From  equations  (4.2),  (4.3),  (4.7)  and  (D. 38)  in  Appendix  D, 
we  conclude  that 


E v(s)  = 0 (4.8) 

and 

E v(si)v(s2)  = 6(s2  - Sj)  2 Qb  Tb  (4.9) 

Comparing  equations  (4.6),  (4.8)  and  (4.9)  with  equations  (D. 1) 
through  (D. 3)  in  Appendix  D,  we  conclude  with  the  aid  of  equations 
(D.  6)  through  (D. 12)  and  (4.4)  and  (4.5)  that 


b(t)  = b ( t o ) + (l-gb)b’(to)  + h(t), 

(4.10) 

where 

b'(t)  * gb  b’(to)  + h’ (t), 

(4.11) 

gb  = exp[(t0-t)/Tp]. 

(4.12) 

and  h(t) , 

h'(t)  are  random  variables  satisfying 

E h(t)  = E h’  (t)  = 0 

(4.13) 

E h(t)2  = Qb[2(t-t0)  - Tb (1-gb)  (3-gb)] 

(4.14) 

E h' (t) 2 - Qb  Tb(l-gb2) 

(4.15) 

E h(t)h' (t)  » Qb  Tb(l-gb) 2 
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(4.16) 

X 


rt 

From  equations  (A. 15),  (4.11)  and  (4.4)  we  deduce  that  for  large 
values  of  t 

E b(t)  2 - Qb/Tb  (4.17) 

From  equations  (D.33),  (D.  34) , (D.  35) , (D.3),  (4.4),  (4.5)  and 
(4.9)  we  find  that  for  large  valuer  of  t 

E_bXt)bXt+4tl  = exp  (4.18) 

E b(t)b(t) 

E[b(t)-b(t+At)] 2 = 2Q[At-Tb(l-exp-At/Tb)],  (4.19) 

If  At  is  small  compared  with  Tb  the  above  equation  reduces  to 

E[b(t)  - b(t+At)] 2 = (Q/Tb) (At) 2 (4.20) 


4. 1 Simulation  of  Clock  Bias. 

Equations  (4.10)  through  (4.16)  are  used  in  the 
simulation  of  clock  bias.  First  h(t)  is  chosen  as  a normally 
distributed  random  number  satisfying  equations  (4.13)  and  (4.14). 
Defining  a,  e,  c by 

a = E h(t) 2 , e = E h’(t)2,  c = E h(t)h’(t),  (4.21) 

h' (t)  is  then  computed  as 

h’(t)  = h(t)  c/a  + k(t),  (4.22) 


-25- 


where  k(t)  is  another  normally  distributed  random  number 
satisfying 


f 


E k(t)  = 0 

and 

E k(t) 2 = e - c 2 /a 


(4.23) 

(4.24) 


Since  h(t)  and  k(t)  are  independently  chosen  random 
variables,  it  follows  that  h(t)  and  h' (t) , computed  as  described 
above,  will  be  consistent  with  equations  (4.13)  through  (4.16). 


-26- 


W 


The  State  Transition  Matrix  and  the  Linearization  of  the 
Equations  of  Motion 


5. 1 Definition  of  terms  and  linearization  of  the  equations 

of  motion. 

Y six  parameter  nominal  position-velocity  vector 

P nominal  (and  constant)  drag  coefficient 

Y + y perturbed  six  parameter  position-velocity  vector 

P + p perturbed  drag  coefficient 

b clock  bias 

b’  clock  bias  rate 

x Nine  parameter  state  vector 


xT  = (yT.  p,  b,  b*)  (5.1.1) 

Y = F(Y,P),  (5.1.2) 

Differential  equation  satisfied  by  nominal 
position-velocity  vector 

Y + y = F(Y  + y,  P + p) , _ (5.1.3) 

Differential  equation  satisfied  by  perturbed 
position  velocity  vector 

Y + y = F(Y,P)  + Fy(Y,P)y  + Fp(Y,P)p,  (5.1.4) 

Linearized  form  of  equation  (5.1.3). 

Fy(Y,P)  is  a 6 x 6 matrix  of  partial  derivatives 
and  Fp(Y,P)  is  a 6-vector  of  partial  derivatives 

<(i(t,s)  = 3y(t)/9y(s),  6x6  state-transition  matrix  (5.1.5) 

D(t , s)  = 3y(t)/8P  (5.1.6) 


where  D(s,s)  = 0.  Vector  of  partial  derivatives 
of  the  6 parameter  state  vector  with  respect  to 
a constant  change  in  the  drag  coefficient  at  time 
s. 

M(t,s)  = 3x(t)/3x(s),  9x9  state-transition  matrix  (5.1.7) 

G(t)  = Fy(Y,P),  H(t)  * Fp(Y,P), 

where  Y is  a function  of  time  but  P is  a constant  (5.1.8) 

e Epoch  Time  (start  time  of  integration  of  nominal  orbit) 
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5.2 

Computational  equations 

for  the 

state  transition 

matrix. 

«• 

* 

M(t,s)  = 

<J>(t,s) 

D(t , s) 

0 

0 

(5.2.1) 

0T 

gp 

o 

o 

0T 

o 

1 

1-gb 

0T 

o 

o 

gb 

. 

where  0 is 

a six 

dimensional  null 

vector 

and  o is 

a scalar 

zero. 

<j»(t,s)  = 4>(t,e)<p(s,e)~  1 (5.2.2) 

D(t,s)  = D(t,e)  - <{>(t , s)D(s ,e)  (5.2.3) 

5.3  Derivation  of  equations. 

The  right-hand  side  of  equation  (5.2.1)  follows  from 
the  definitions  of  M(t,s),  4>(t,s)  and  D(t,s)  and  also  from 
equations  (3.10),  (4.10)  and  (4.11).  Equations  (5.2.2)  and 

(5.2.3)  remain  to  be  derived.  It  follows  from  equations  (5.1.2), 

(5.1.4)  and  (5.1.8)  that 

y(t)  - G(t)y (t)  + H(t)p (t) . (5.3.1) 

<J>(t,e)  is  obtained  as  the  solution  of  the  differential  equation 

4>(t, e)  = G(t)4»(t,e)  with  4>(e,e)  = I,  (5.3.2) 

and  D(t,e)  as  the  solution  of 
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D(t,e)  = G(t)D(t , e)  + H(t), 


(5.3.3) 


with  D(e,e)  = 0 


We  wish  to  show  that 


(5.3.4) 


y(t)  = <j>(t,  s)y(s)  + D(t,  s)p(s)  (5.3.5) 

is  a solution  of  equation  (5.3.1)  and  furthermore  that  in  order 
that  the  left  and  right  hand  sides  of  equation  (5.3.5)  be 
consistent. 


4* (s , s)  = I and  D(s,s)  = 0 (5.3.6) 

Differentiating  equation  (5.3.5)  with  respect  to  t we  deduce 
with  the  aid  of  equations  (5.2.2),  (5.2.3),  (5.3.2)  and  (5.3.3) 
that 

y (t)  = G(t)<{>(t,  s)y(s)  + [ G(t)D(t , e)  + H(t) 

- G(t)0(t,s)D(s,e)]p(s)  , 

i.  e. , 

y (t)  = G(t)<j>(t,  s)y(s)  + [ G(t)D(t , s)  +H(t)]p(s) 

- G(t)[<Kt,s)y(s)  + D(t,s)p(s)]+  H(t)p(s)  . 

With  the  aid  of  equation  (5.3.5)  the  above  equation 
can  be  seen  to  reduce  to 
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y(t)  = G(t)y(t)  + H(t)p(s) 


(5.3.7) 


Equations  (5.3.6)  can  be  seen  to  follow  from  equations 
(5.2.2)  and  (5.2.3).  Equation  (5.3.7),  however,  is  not  identical 
to  equation  (5.3.1),  but  it  is  a good  approximation  to  it  pro- 
vided that  the  drag  perturbation  p changes  but  little  in  the  time 
interval  (s,t).  This  we  shall  assume  to  be  the  case. 


6.  Simulation  of  Measurements  Involving  the  Global  Positioning 
System  (g!PS) 

6. 1 Brief  Description  of  GPS 

GPS  consists  of  a set  of  satellites,  whose  positions 
and  velocities  are  known  to  all  users  of  the  system.  These 
satellites  transmit  radio  signals  at  fixed  intervals.  The 
clocks  of  the  GPS  satellites  are  extremely  accurate.  They  are 
also  mutually  synchronized.  If  the  user's  clock  also  were 
synchronized  with  the  GPS  clocks,  then  the  user  could  calculate 
his  distance  to  each  GPS  satellite  (provided  of  course  that 
he  could  see  it).  Given  three  distances  to  three  known  positions, 
the  user  may  then  solve  a simple  geometric  problem  to  obtain  his 
own  position.  If  the  user  clock  is  not  very  accurate,  then  the 
user  may  instead  process  the  signal  from  a fourth  satellite  to 
give  similar  results. 

Specifically,  GPS  consists  of  24  satellites  arranged 
in  3 rings  of  8 equally  spaced  satellites  (see  Figures  6.1  and 
6.2).  Each  satellite  is  in  a 12  hour  (26610  km  radius)  circular 
orbit  with  an  orbital  inclination  of  63  degrees.  The  longitudes 
of  the  ascending  nodes  of  the  satellite  orbits  are  0 degrees 
for  those  in  ring  1,  120  degrees  for  those  in  ring  2,  and  240 
degrees  for  those  in  ring  3.  Since  the  satellites  of  each  ring 
are  equally  spaced  the  angular  distance  between  them  must  be 
45  degrees.  For  each  ring  a satellite  must  thus  cross  the 
equator  from  South  to  North  every  90  minutes  (another  satellite 
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simultaneously  crosses  from  North  to  South).  The  satellites  of 
the  three  rings  are  phased  relative  to  each  other  such  that  a 
satellite  will  cross  the  equator  South  to  North  every  30  minutes. 

The  order  is  ring  1,  ring  2,  ring  3,  ring  1,....  The  three  rings 
must  obviously  intersect  one  another..  However,  no  two  satellites 
will  ever  approach  each  other  closer  than  10.4  degrees.  The 
orbital  paths  intersect  each  other  at  a latitude  of  44.5  degrees 
(North  and  South).  At  the  point  of  intersection,  satellites  of 
two  different  rings  approach  each  other  at  101  degrees.  The 
longitudes  of  the  intersections  in  the  Northern  Hemisphere  occur 
at  30  degrees  (1  ascending,  3 descending),  150  degrees  (2  ascending, 
1 descending)  and  270  degrees  (3  ascending,  2 descending). 


6.2  Cartesian  Coordinates  of  the  GPS  Satellites 

The  position  of  each  GPS  satellite  may  be  specified 
(i)  by  the  longitude  ft  of  the  ascending  node  of  its  orbit 
(0  degrees  for  ring  1,  120  degrees  for  ring  2,  240  degrees  for 
ring  3),  and  (ii)  by  its  angular  distance  to  from  that  node,  w 
is  computed  from  the  formula 


to  = toQ  + nAto  + <ot, 


(6.2.1) 


where  <oQ  equals  0 degrees  for  ring  1,  30  degrees  for  ring  2,  and 

15  degrees  for  ring  3, 

A to  equals  45  degrees , 

n is  the  satellite  number 

(n  = 0 , 1 , ....  7 for  each  ring) , 

to  = 360  degrees/12  hours,  and  t is 
time  from  midnight. 


i 


A 
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Let  p and  q be  two  unit  vectors  lying  in  the  orbital  plane, 
p pointing  towards  the  ascending  node  and  q pointing  towards 
the  point  of  highest  latitude.  Then 

p = (cos  SI,  sin  Si,  0)  and  (6.2.2) 

q = (-sin  SI  cosi,  cos  ft  cos  i,  sin  i)  , (6.2.3) 

where  i is  the  orbital  inclination.  The  satellite  position  r 
is  then  given  by 

r = p cos  ui  + q sin  w (6.2.4) 

6.3  GPS  Measurements  and  their  Partial  Derivatives 

The  Cartesian  coordinates  of  each  GPS  satellite  is 
given  by  a formula  of  the  form  (6.2.4).  To  distinguish  between 
the  different  satellites  we  add  a subscript.  Thus,  r^ , 

(j  = 1,  2,  ...24)  is  defined  as  the  position  vector  of  the  j-th 
GPS  satellite.  Similarly  we  define  r as  the  position  vector  of 
the  user.  The  GPS  measurement  to  GPS  satellite  number  j is  then 
given  by 

dj  = ^3  " + b*  (6.3.1) 

where  b is  a bia3  term  due  to  a user  clock  error.  We  define  the 
unit  vector  v.  by 

Vj  = (£j  " t)/(dj  " b)  (6.3.2) 
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Differentiating  with  respect  to  r we  obtain 


3dj/8r  = -Vj 

Also 

3d.. /3b  = 1 


(6.3.4) 


The  GPS  measurement  vector  (Z(k)  in  Section  2)  is 
made  up  of  four  measurements  of  the  form  (6.3.1).  The  partial 
derivative  matrix  (H(k)  in  Section  2)  is  made  up  of  the  corres- 
ponding partial  derivatives  as  given  by  equations  (6.3.3)  and 
(6.3.4). 


6.4  GPS  Simulations 

The  problem  to  be  solved  in  GPS  simulations,  just  as 
in  real  situation  scheduling,  is  how  to  choose  4 GPS  satellites 
out  of  24  so  as  to  be  able  to  derive  the  best  possible  user 
position.  Since  a user  satellite  has  a clear  'horizon',  he  can 
see  roughly  a hemisphere  of  GPS  satellites.  This,  on  the  average, 
amounts  to  12  satellites.  There  are  495  different  ways  to  pick 
4 out  of  12.  To  find  the  best  4 it  is  necessary  to  test  each 
combination.  To  do  so  at  every  time  point  is  impractical.  The 
following  is  a suboptional  but  good  scheme  for  the  selection 
process.  (Further  details  are  given  in  Appendix  E.) 

6.4.1.  Satellite  visibility.  It  is  of  course  necessary 
that  each  selected  GPS  satellite  be  visible  to  the  user.  In  order 
that  the  satellite  not  be  visible  two  criteria  must  be  met 
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(i)  The  satellite  must  appear  below  the  user’s 
'horizon' , i.e. 


VjT  r < 0, 


(6.4.1) 


(ii)  The  satellite  to  user  line  of  sight  must 

intersect  the  Earth,  defined  for  this  purpose 
as  including  an  atomsphere  100  km  above  the 
surface,  i.e. 


rTr  < (VjTr ) + re2 


(6.4.2) 


where  r is  the  radius  of  the  earth  as  defined 
above . 

6.4.2.  Selection  of  first  satellite.  The  first 
satellite  is  somewhat  arbitrarily  chosen  as  the  one  that  is 

T_ 

highest  in  the  sky  (v^  r is  a maximum) . There  is  no  loss  of 
generality  in  designating  this  as  satellite  number  one  (j  = 1) . 

6.4.3  Selection  of  second  satellite.  It  can  be 
shown  (see  Appendix  E)  that  the  optimal  geometric  configuration 
obtains  when  the  angles  between  the  four  lines  of  sight  are  all 
equal.  This  is  possible  only  if  the  angles  equal  cos-1 (-1/3)  or 
109.5  degrees.  The  second  satellite  is  therefore  chosen  such 
that 

IviS  + 

is  a minimum.  This  satellite  is  designated  satellite  number  2 

(j  - 2). 

6.4.4  Selection  of  the  third  satellite.  It  is  shown 

in  Appendix  E that  given  two  satellites  (1  and  2)  then  the  optimal 
geometric  lines  of  sight  for  satellites  3 and  4 must  satisfy: 
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(i)  v3  and  v4  lie  in  a plane  perpendicular 
to  that  defined  by  Vi  and  v2 , 

(ii)  v3  + v4  is  diametrically  opposite  to  Vi  + v2. 

(iii)  If  the  angle  between  v3  and  v2  is  2a,  and  the 
angle  between  v3  and  v4  is  2$,  then 


cos3  = .327  cosa  - .765 


(6.4.2) 


The  procedure  then  is  to  define  two  unit  vectors  u3  and  u4 

satisfying  the  above  three  criteria  and  then  finding  the  third 

T 

satellite  such  that  vju*  where  u = u3  or  u4,  is  maximized. 

This  satellite  is  designated  number  3 (j  = 3) . u3  and  u4  are 
computed  as  follows: 

cos  a = fix  + v?  v2)/2  (6.4.3) 


cos3  is  then  computed  using  equation  (6.4.2). 
Hence 

sin  3 = /I  - cos 2 3 


The  vector  c is  defined  by 

c = (Vi  + v2)  cos3/2cosa, 
and  the  vector  e by 

e - (vi  a v2)sin3[l  - (v?vz)2J  * 


Then 


u3  «»  c + e and  u4  = c - e 


(6.4.4) 


(6.4.5) 


(6.4.6) 


(6.4.7) 


1 


-36- 


6.4.5  Selection  of  the  fourth  satellite.  Given  3 
satellites  the  fourth  one  is  selected  optimally  as  follows: 
(For  derivation  of  the  formulae  see  Appendix  E).  Define  the 


matrix 


Then 


T 

A0  = [vi,  v2,  v3]  . 


(6.4.8) 


Ao1  = [v2av3,  v 3 av i , ViAv2]/ [(vjAv2)*v3]  (6.4.9) 


Define  the  vector  d0  by 


do  - [1,  1,  1]^ 


Then  compute 


(6.4.10) 


T 

e0  = Ao  1 dD,  and  e00  = Ao  eo 


(6.4.11) 


The  fourth  satellite  is  then  found  by  minimizing  the  expression 


where 


2 hi,  e00T  f + hi  e0Te0  (l+fTf ) , 


v-  , 


(6.4.12) 


h4  = (1  - v4e0) 


(6.4.13) 


f = AoTv„ 


(6.4.14) 
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GPS  SATELLITE  SYSTEM 

AS  VIEWED  FROM  A DISTANT  POINT  AT  30°  LATITUDE 


FIGURE  6.2 
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GPS  SATELLITE  SYSTEM 

AS  VIEWED  FROM  A DISTANT  POINT  ABOVE  THE  NORTH  POLE 
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7.  Drag  Segments  in  an  Orbital  Arc. 

In  most  cases  the  drag  coefficient  of  an  orbital 
satellite  is  a constant.  This  is  true  even  if  the  shape  of  the 
satellite  is  not  spherically  symmetric,  provided  that  the  satellite 
presents  the  same  aspect  angle  along  its  direction  of  motion. 

In  the  equations  of  motion  the  drag  coefficient  always  appears 
as  a factor  multiplying  the  atmospheric  density.  It  is  thus 
possible  to  compensate  for  density  variations  through  corresponding 
changes  in  the  drag  coefficient.  This  is  often  done  in  practice. 
The  capability  to  do  that  has  now  been  added  to  Photonap.  The 
program  has  been  modified  to  include  a number  of  different  drag 
segments.  Each  drag  coefficient  may  be  constrained  either  to 
some  a priori  value  (absolute  constraint)  or  the  coefficients 
of  contiguous  segments  may  be  constrained  relative  to  each  other 
(relative  constraints). 


7 . 1 Partial  Derivatives  of  the  State  Vector  with  Respect 

to  the  Drag  Coefficients. 

Let 

Y=F(Y,Pk),  (7.1) 

denote  the  differential  equation  governing  the  satellite  motion. 

Y is  the  state  vector  (position-velocity  vector)  and  Pk  is  the 
drag  coefficient  in  segment  k valid  in  the  interval  between  times 
tk  and  tk+p  Let  y and  pk  denote  a small  perturbations  in  Y and 
Pk,  respectively.  Then 
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y - F(Y.Pk)  y + F(Y,Pk)pk 

(7.2) 

Writing 

G « 3F(Y,Pk)/3Y 

(7.3) 

and 

H = 3F(Y,Pk)/3Pk 

(7.3) 

the  above  equation  may  be  written  as 


y - Gy  + Hpk  (7. A) 

To  indicate  that  y,  G and  H are  functions  of  time  we 
rewrite  the  equation  in  the  form 

y (t)  - G(t)y (t)  + H(t)pk  (7.5) 

The  6x6  state  transition  matrix  <f>(t,e)  is  defined 
as  the  solution  of  the  differential  equation 

<j>(t,e)  = G(t)«Kt,e),  (7.6) 

with 

<Ke,e)  = I,  (7.7) 

e being  the  time  of  Epoch. 

The  six-vector  of  drag  partials,  Dk(t)  is  defined  as 
the  solution  of  the  differential  equation 

Dk(C)  = G(t)Dk(t)  + H(t),  (7.8) 


with  D.  (t,  ) = 0, 


(7.9 


and  Dj^t)  being  defined  only  in  the  interval  [t^,  tk+1j.  Equation 
(7.5)  is  similarly,  for  each  value  of  k,  valid  only  in  the  interval 
[t^,  tic+l]'  ^e  so^uti°ns  °f  equations  (7.5)  must  be  continuous 
and  satisfy  the  initial  conditions 

y(ti)  = o,  (7.10) 

where  ti  = e,  the  Epoch  Time. 

It  will  now  be  shown  that  if 

Ck  5 C S ck+l  {7.11) 

then  equation  (7.5)  is  satisfied  by 

k 

xk(t)  = <p(t,e)  <J>( ti,e)-1Di_1(ti)  p±_1 

+ Dk(t)pk  (7.12) 

Differentiating  equation  (7.12)  with  respect  to  t, it 
can  easily  be  seen  that  equation  (7.5)  is  satisfied.  It  remains 
to  be  shown  that  the  solution  is  continuous,  i.e.,  that  *k(tk)  = 
xk_1(tk).  From  equation  (7.12)  we  obtain 

k-1 

*k-l(tk>  * ♦<tk-e>  «ti •e>“ 

+ Wi^k-l-  <713> 

The  above  equation  may  be  rewritten  as 


k 

xlc-1  = 4>(ti.,e)  S <(>(t.  ,e)-1  D.  ,(t.)p. 


(7.14) 


From  equations  (7,9),  (7.12)  and  (7.14)  we  deduce 


that 


xk-i(tk>  ■ xk(tk> 


(7.15) 


Since  xi(ti)  = 0 it  follows  that  x^(t)  as  given  by  equation 
(7.12)  is  the  required  solution  of  equation  (7.5),  i.e.,  for  t 
satisfying  equation  (7.11), 


y(t)  = Kt,e)  2 Kt^.e)-1 

+ Dk<c)pk 


(7.16) 


It  hence  follows  that  the  required  partial  derivatives 


are  given  by 


for  t > ti+1,  8y(t)/3pi  = *(t,e)Mti+1,e)"1D.(ti+1) 


for  < t < ti+1,  3y(t)/3pi  = Di(t) 


(7.17) 


for  t < t^,  3y(t)/3pi  = 0 


7.2  Absolute  and  Relative  Constraints 

To  obtain  the  solution  6 p of  a linearized  weighted  least 
squares  problem,  an  equation  of  the  following  form  must  be  solved 


N 6p  - b 


(7.18) 
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In  the  above  equation,  N is  a positive-definite  matrix 


(the  normal  equations  coefficient  matrix)  and  b a vector  (the 
normal  equations  vector) . Using  the  summation  convention  equation 
(7.18)  may  be  written  in  index  form  as 

N(i.j)6p(j)  = b(i)  (7.19) 

If  all  measurements  are  statistically  independent,  then 
N(i,j)  and  b(i)  are  computed  as  the  sum  of  terms  of  the  form 

(7.20) 

and 

Ab(i)  = W (mo-m).  (7.21) 

where  W,  the  measurement  weight  is  inversely  proportional  to  the 
variance  of  the  measurement  error, 

m0  is  the  observed  measurement 
and  m is  the  measurement  calculated  as  a function 
of  a set  of  parameters  p(k). 

After  the  solution  of  equation  (7.19)  has  been  obtained 
the  estimated  parameter  is  updated  to  p(k)  + 5p(k). 

7.2.1  Absolute  Constraints.  If,  a priori,  we  know  that 
parameter  p(i)  = a^  and  that  the  error  in  a ^ has  a variance  of  a2, 
then  we  may  treat  that  information  in  exactly  the  same  way  as 
measurement  information.  Hence  mo  = a^,  m ® p(i)  and  3m/8p(i)  = 1. 


-44- 


In  accordance  with  equations  (7.20)  and  (7.21)  the  contributions 
to  N and  b are  then  given  by 


AN(i,i)  = o'2  (7.22) 

Ab(i)  = o"2(ai  - p(i) ) (7.23) 

7.2.2  Relative  Constraints.  If,  a priori,  we  know 
that  parameters  i and  j should  assume  the  same  value  and  that  the 
error  in  this  assumption  has  a variance  of  a2,  then  we  may  treat 
this  information  exactly  the  same  way  as  measurement  information. 
Hence  m0  = 0,  m = p(i)  - p(j),  9m/9p(i)  = 1 and  9m/9p(j)  = -1. 

In  accordance  with  equations  (7.20)  and  (7.21)  the 
contributions  to  N and  b are  then  given  by 


and 


AN(i,i)  = a" 2 
AN(i,j)  = -a- 2 
AN(j  ,j)  = a~ 2 


Ab(i)  = a" 2 p(j)  - p(i) 
Ab(j ) = a"2  p(i)  - p(j) 


(7.24) 


(7.25) 


Note  that  m could  equally  well  have  been  defined  by  m = p(j)  - p(i) 
The  result,  however,  would  be  the  same. 
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The  Lockheed- Jacchia  Atmospheric  Model 


The  equations  presented  in  this  section  are  based 
partly  on  some  equations  supplied  by  Mr.  George  Stentz  of  DMAAC, 
St.  Louis,  MO  and  partly  on  a computer  program  listing  from  the 
same  source.  Part  of  the  description  comes  from  (Jacchia,  1960). 
A modified  version  of  the  DMAAC  program  has  been  incorporated  in 
Photonap . 


Description  of  Variables  and  Constants  used  in  the  Program. 


Angle  between  point  of  interest  and  point  of  maximum 
solar  heating  effect  as  seen  from  the  center  of  the 
Earth. 

_ (l+cosv \ ■ 

= .55  radians.  Lag  angle.  Angle  between  the  sun  and 
the  point  of  maximum  solar  heating  as  seen  from  the 
center  of  the  Earth. 

Time  (in  days)  from  noon  on  January  1,  4713  B.C. 

= 2436204.  Number  of  days  between  January  1,  1958 
and  January  1,  4713  B.C. 

Time  (in  days)  since  noon  on  December  31,  1957. 

= .017203  radians/day.  The  Earth’s  orbital  rate 
about  the  sun 

Longitude  of  the  Sun  measured  along  the  ecliptic  from 
the  equinox 

-Longitude  of  the  Sun  when  t ' = 0 
= 0.01675.  The  eccentricity  of  the  Earth’s  orbit 


,4092  radians.  Obliquity  of  the  ecliptic 
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1 n.m. 


10.7  cm  flux  mea  ured  in  units  of  100  x 10~22 
watt/m2/ cycle/ sec 

20  cm  flux  measured  in  same  units  as  Fj0.7 

_ Q*gj  constants  used  in  computation  of  F10.7 
2 71 

= T7J27)  Frequency  corresponding  to  a period  of 
4020  days  (approximately  11  years) 

= 0.85  Conversion  factor  for  converting  F10.7 
to  the  equivalent  F20 

Atmospheric  Density  (slugs/cu. f t) 

" p = HE  log  p 

Height  above  the  surface 

= 76  n.m. 

= 108  n.m. 

= 378  n.m. 

*=  1000  n.m. 

= 5.606  x 10" 1 2 slugs/cu. ft. 

= 7.18  unless  otherwise  specified  by  user 
= 153  n.m. 

= -15.738  unless  otherwise  specified  by  the  user 
= .00368  (n.m.)-1 

- 6.363 

= .0048  (n.m.)-1 
= 0.19 

= .0102  (n.m.)’1 

- 1.9 

= .00504  (slugs/cu. ft. ) (n.m. ) 5 
* 6 x 106  (n.m. ) 3 

“ 1.852  km 


1 slug/cu.ft.  = 0.515378  x 1012  kg/knr 
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8. 2 Description  of  Equations 

The  Jacchia  Atmospheric  Model  is  a dynamic  model  in 
the  sense  that  it  is  a function,  not  only  of  position,  but  also 
of  time.  The  time  dependence  is  due  to  solar  heating.  However, 
since  solar  heating  is  not  instantaneous,  the  maximum  perturba- 
tion to  the  atmosphere  will  occur  some  time  after  noon,  local 
time  (according  to  this  model  just  after  2 p.m.  local  time).  If 
5 is  a unit  vector  pointing  towards  the  sun,  then  the  maximum 
perturbation  will  occur  in  the  direction  of  s',  where  s and  s' 
point  towards  the  same  latitude,  but  s'  towards  a point  X radians 
further  East.  If  u is  a unit  vector  pointing  towards  the  point 
of  interest,  then  cos  ip  is  defined  by 


cos  tji  = u • s' 


If 

u = (x,y , z) 

and 

S = (si,s2,s3), 

then 

s'  = (sicosA  - s2sinX,  s2cosX  + SisinX,s3) 


and 

costy  = (six  + s2y)cosX  - (s2x  - siy)  sinX  + s3 
5 is  computed  using  the  following  equations 


t* 


L = wt*  + 2e  sin  u>t’  - Lo 
5 = (cosL,  sinL  cose,  sinL  sine) 


(8.1) 

(8.2) 

(8.3) 

(8.4) 

(8.5) 

(8.6) 

(8.7) 

(8.8) 
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Combining  equations  (8.3),  (8.5)  and  (8.8)  we  obtain 


cosip  = (xcosL  + ysinL  cosc)cosA 

+ (-xsinL  cost.  + ycosL)sinA  + sinLsine  (8.9) 

Let  g be  defined  by 

g = [cosi|j/2]6  , 

i.  e.  , 

g . (l+cosl)*  (8.10) 

Unless  input  by  the  user  the  10.7  cm  flux  is  given  by 

K i o . 7 = Cg  + C9  + cos(u)pt')  (8.11) 

This  is  converted  to  an  equivalent  flux  at  20  cm  by  the  formulae 

F20  = C15  K 1 0 . 7 (8.12) 

For  the  purpose  of  calculating  the  density,  the  atmosphere  is 
subdivided  into  four  different  regions,  with  different  sets  of 
formulae  being  valid  in  each  region.  These  are  given  below. 


Region  A 

hi < h < h2 

P - P 1 P2 P 3 » 

(8.13A) 

where 

p.  - -&)d‘ 

(8.14A) 

pi  ’[kFK,  + F”] 

(8.15A) 
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P 3 


1 + lizlli 

Cl  8 


g 


(8.16A) 


Region  B 


where 


and 


• _ di  . 1 I"  1 i ^Fio  T.l 

p "T+pT  '[iTFKT  + E7^K7j  + 

h2<  h < h3 
p = Po  q 

Po  = 10dz  “ C2eh  + °27  exP(_c2eh) 


S . 

Ci  aPa 


P* 


q = F20|l  + c29  £exp(c30h)-c3  ijgj 
loge10[-c26-c28c27exp(-c28h)]+  .F2  oC29C?^,xp(c30h)g 


(8. 17A) 


(8.13B) 

(8.14B) 

(8.15B) 


Region  C h3<h  < h* 


where 


p = b ib2 , 


bi 

b2 


- C 3 6 


Fio. 


gd  - fF  ) + 


C 3 7 


5 3/7  _\  c 3 


- E - K 


h’  b. 


Region  D h i hi,  or  h < hi 
P - p’  - 0 


(8.13C) 

(8.14C) 

(8.15C) 

(8. 16C) 

(8.13D) 


Note  that  in  all  of  the  above  equations  the  density  is  computed 
in  slugs/cu.ft.  and  p'  in  1/n.m.  Before  being  used  by  Photonap 
these  quantities  are  converted  to  kg/km*  and  km”1,  respectively 
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9.  An  Outline  of  Program  Changes 

The  program  changes  made  to  Photonap  fall  into  two 
categories:  (i)  changes  to  existing  routines,  and  (ii)  the 

addition  of  new  routines.  Two  of  the  existing  routines,  SPOLCD 
and  SOLVER,  were  initially  simply  modified,  but  owing  to  the  routines 
in  the  process  becoming  extremely  lengthy  and  unmanageable,  they 
were  later  split  into  smaller  routines.  Thus  SPOLCD  was  split 
into  SPOLCD,  SP00L1,  SP00L2  and  SP0100.  SOLVER  was  split  into 
SOLVER,  S0LV1,  S0LV2  and  S0LV3.  The  following  totally  new 
routines  have  been  added: 

A.  Routines  associated  with  the  Lockheed- Jacchia 

atmospheric  model:  JACHIA 

B.  Routines  associated  with  drag  segmentation:  DRAGU 

C.  Routines  associated  with  Kalman  filtering  and  smoothing: 

KMNCON,  KMNEVA,  KMNIDE,  KMNINI,  KMNINV , KMNMP1,  KMNMP2, 
KMNMP3,  KMNRAN,  KMNOUT,  KMNSIM,  KMNSM2 , INVSYM,  INVSYS, 
MAT99,  VXPROD . 

D.  Routines  associated  with  GPS  measurements:  SELECT 

E.  Routines  associated  with  normal  equations  for  correlated 
measurements:  SOLVFU,  SOLAWA 

In  addition  to  the  changes  described  above,  subroutine 
SVARED,  after  a trivial  change  in  the  coding,  was  found  to  be 
superfluous,  and  was  hence  removed  from  Photonap. 

A flow  chart  of  Photonap  together  with  a short  descrip- 
tion of  each  routine  is  given  in  Appendix  F. 


A 
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10. 


Changes  to  the  Photonap  User's  Guide 


(i) 


(ii) 

(iii) 

(iv) 


(v) 

(vi) 


(vii) 


(viii) 


(ix) 


The  following  has  been  added  to  the  User's  Guide. 
Insertion  into  Section  I of  a general  description 
of  control  card  set-ups  for  Kalman  filtering  and 
smoothing. 

Addition  to  101  card  input  to  specify  Kalman  filter  mode. 
Description  of  230  card  for  specifying  drag  segments. 
Addition  of  note  (Note  13)  to  601  card  for  handling  of 
drag  coefficients  appearing  in  different  drag  segments 
of  the  same  arc. 

Description  of  612  card  for  specifying  constants  required 
by  Lockheed- Jacchia  atmosphere 

Description  of  614  card  for  specifying  GPS  filter 
constants 

Addition  of  note  (Note  5)  to  701  card  for  processing  of 
correlated  measurements  output  from  Kalman  filter  or 
smoother. 

Addition  to  Appendix  IB  describing  tape  format  for 
Kalman  filter  input,  and  tape  format  for  Kalman  filter 
output . 

Addition  of  Appendix  V describing  example  of  the  job 
control  language  required  for  running  Photonap  on  the 
CDC  6400. 
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11.  Test  Runs 

In  order  to  check  out  the  modified  version  of  Photonap 
a large  number  of  test  runs  were  made.  These  included  running 
the  standard  set  of  Photonap  test  decks,  which  are  run  after 
every  program  modification.  Five  new  test  decks,  designated  PB4, 
PB5,  PB6,  PB7  and  PB8,  have  been  added  to  the  standard  set,  which 
now  consists  of 

non-photogrammetric  test  decks  TESTXX,  PA1,  PA2, 

PA3,  PA4,  PA5,  PA5X1,  PA5X2 , PA6 , PA7,  PA8,  PA9 , 

PBO,  PB1,  PB2,  PB3,  PB4,  PB5 , PB6,  PB7,  PB8, 

Photogrammetric  test  decks  PAA,  PAB,  PAC,  PAD, 

PAE,  PAF,  PAG,  PAGX1 , 

combined  test  deck  FATBOY. 

A short  description  of  each  of  the  new  test  decks  is  given  below. 

11.1  Test  Deck  PB4.  Lockheed- Jacchia  atmosphere  and  multiple 

drag  segments.  Two  parts. 

(i)  Data  generation  using  Jacchia  Atmosphere  and  a 
single  drag  segment. 

(ii)  Orbit  and  drag  coefficient  recovery  using  U.S. 
Standard  Atmosphere.  Six  drag  segments  with 
relative  constraints.  Epoch  coincident  with  start 
of  first  drag  segment. 
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11.2  Test  Deck  PB5.  Lockheed- Jacchia  Atmosphere  and  multiple 
drag  segments.  Two  parts. 

(i)  Data  generation  using  Jacchia  Atmosphere  and  a 
single  drag  segment. 

(ii)  Orbit  and  drag  coefficient  recovery  using  U.S. 

Standard  Atmosphere.  Six  drag  segments  with 
relative  and  absolute  constraints.  Epoch  in 
middle  of  fourth  segment. 

11.3  Test  Deck  PB6.  Six  point  smoother  using  GPS  measurements. 
Three  part  run. 

(i)  GPS  data  generation  using  Lockheed- Jacchia  Atmosphere, 

(ii)  Six  point  smoother  using  U.S.  Standard  Atmosphere. 
Smoother  output  of  position  and  velocity  at  30 
second  intervals. 

(iii)  Orbit  comparison  between  the  smoother  output  and 
the  orbit  used  in  data  generation. 

11.4  Test  Deck  PB7.  Filter  using  GPS  measurements.  Three 
part  run. 

(i)  GPS  data  generation  using  Lockheed-Jacchia  Atmosphere, 
(ii)  Filter  (0  point  smoother)  using  U.S.  Standard 

Atmosphere.  Filter  output  of  position  and  velocity 
at  30  second  intervals. 

(iii)  Orbit  recovery  based  on  Lockheed-Jacchia  Atmosphere. 
Filter  output  used  as  measurement  data. 
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11.5  Test  Deck  PB8.  Filter  and  smoother  comparisons. 

U.S.  Standard  Atmosphere  used  in  all  three  parts: 

(i)  GPS  data  generation, 

(iia)  Filter, 

(iib)  16  point  smoother. 


Comparison  between  orbit  used  in  generation 
(Part  (i) , pages  7 through  9)  and  the  orbit 
recovered  (Parts  (ii) , pages  4 through  6) 
shows  the  superiority  of  the  16  point  smoother 
over  the  filter. 
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APPENDIX  A 


Some  formulae  for  the  differentiation  of  the  trace  of  matrix 


products . 

Given  a positive-definite  matrix  D and  general  matrices 
A and  F,  the  following  three  formulae  will  be  derived. 


TF  Trace 

(DFAFT)  = DF(A  + AT) 

(A.l) 

•gy  Trace 

(DFA)  = DAT 

(A.  2) 

Trace 

(DAFT)  - DA 

(A.  3) 

Proof . 


Denoting  the  left  hand  side  of  formula  (A. 1)  by  L, 
we  may  express  it  in  index  form  as 


Lij  3Ftj 


£ 

a,b,c,d 


DabFbcAcdFad 


- h IWaa  + & 


= (DTFAT  + DFA)^ 


Since  D is  symmetric,  formula  (A.l)  follows.  Denoting 
the  left  hand  side  of  the  second  formula  by  M,  we  similarly 
obtain: 
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Mij  " E DabFbcA 

J a,b,c 


ca 


- £ D .A. 
a ai  3a 


= (DTAT) . . 

ij 


Since  D is  symmetric,  the  second  formula  follows. 
Denoting  the  left  hand  side  of  the  third  formula  by  N,  we  find 
that 


N 


ij  " WT  E DabAbcF 


ij 


a,b,  c 


ac 


^ °ibAbj 


(DA)ij ’ 


which  is  equivalent  to  the  third  formula. 
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where  I is  the  identity  matrix.  Denoting  the  errors  in  x,  a 
and  b by  x,  a and  b,  respectively,  we  deduce  from  equations 
(B.3)  and  (B.4)  that 

x - Fa  + (I-F)b  (B. 5) 


Since  a and  b are  independent  it  follows  from  the  above  and  the 
definition  of  the  covariance  that 

P = FAFT  + (I-F)B(I-F)T  (B. 6) 

T 

F is  chosen  such  that  the  expected  value  of  x Dx,  where  D is  a 
positive-definite  matrix,  is  minimized.  It  turns  out  that  as 
long  as  D is  symmetric  and  non-singular  the  solution  is  independent 
of  the  choice  of  D.  Remembering  that  if  XY  and  YX  are  both 
square  matrices,  then  trace  (XY)  = trace  (YX) , it  follows  that 
the  quantity  we  are  trying  to  minimize  is  the  expected  value  of 
trace  (Dxx  ),  i.e.,  trace  (DP).  From  equation  (B.6)  and  formulae 
(A.l),  (A. 2)  and  (A. 3)  in  Appendix  A,  we  then  deduce  that 


D[2F(A  + B)  - 2B ] = 0. 

Hence 

F(A  + B)  - B - 0, 

and 

F » (A"  1 + B“  1 )~  'A-  1 


(B.7) 

(B.8) 

(B.9) 


i 
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From  the  above  and  equation  (B.4)  we  find  that 


F1  = (A-1  + B“ *)" XB“ 1 , 


(B. 10) 


and  from  equations  (B.6)  and  (B.8)  we  obtain, 


P = B FB 

- BF1 

From  the  above  and  equation  (B.10)  it  follows  that 
P = (A-  1 + B"  1 ) ~ 1 , 

which  is  equation  (B.2).  Equation  (B. 1)  follows  from  the  above 
and  equations  (B.3),  (B.9)  and  (B.10). 
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APPENDIX  C 


Updating  the  minimum  variance  estimate  based  on  new  independent 
measurements . 


Given 

(i)  an  unbiased  estimate  a of  a parameter  vector, 
whose  true  value  is  xT, 

(ii)  A,  the  covariance  of  the  error  in  a, 

(iii)  a measurement  vector  Z satisfying  the  equation 


Z = H xT  + r,  where  (C.l) 
E (r)  = 0 E(rrT)  = R,  (C.2) 
and  H is  a given  matrix, 

then  the  new  minimum  variance  estimate  x is  given  by 


x = a + G(Z-Ha) , 

where 

G = AHT(R  + HAHT)_1 , 

and  P,  the  covariance  of  the  error  in  x,  is  given  by 


(C.3) 

(C.4) 


P - A - GHA 


(C.5) 


Proof.  Equation  (C.3)  is  clearly  a general  form  for  a linear 
unbiased  estimate  of  x.  We  shall,  therefore,  assume 
that  x is  given  by  equation  (C.3)  and  then  proceed  to 


derive  equations  (C.4)  and  (C.5).  Let  x and  a denote  the  errors 
in  x and  a,  respectively.  It  then  follows  from  equations  (C.l), 
(C.2)  and  (C.3)  that 


x = a + G (r  - Ha),  i.e. 


x = (I  - GH)a  + Gr 


(C.6) 


Since  a and  r,  by  assumption,  are  independent,  it  follows  that 

P = (I-GH)A(I-GH)T  + GRGT  (C. 7) 

G is  chosen  such  that  the  expected  value  of  x Dx,  where  D is  a 
positive- definite  matrix,  is  minimized.  It  turns  out  that  as 
long  as  D is  symmetric  and  non- singular  the  solution  is  independent 
of  the  choice  of  D.  Remembering  that  if  XY  and  YX  are  both  square 
matrices,  then  trace  (XY)  = trace  (YX) , it  follows  that  the  quantity 
we  are  trying  to  minimize  is  the  expected  value  of  trace  (Dxx  ) , 
i.e.,  trace  (DP).  From  equation  (C.7)  and  formulae  (A. 1) , (A. 2) 
and  (A. 3)  in  Appendix  A,  we  then  deduce  that 

D[2G(HAHT+R)  - 2AHT]=  0 


Hence , 


G(HAHT+R)  - aht  = 0 


(C.  8) 


from  which  equation  (C.4)  immediately  follows.  Post-multiplying 

T 

equation  (C.8)  by  G and  subtracting  the  result  of  equation  (C.7) 
yields  equation  (C.5).  This  completes  the  proof. 
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APPENDIX  D 

Differential  equation  associated  with  timing  bias  and  variations 
in  the  drag  coefficient. 

Given  the  differential  equation 

x(t)  + £(t)  = v(t),  (D.l) 

where  a dot  denotes  differentiation  with  respect  to  t,  and  v(t) 
is  a random  variable  with 

(D.2) 
(D.3) 

(D.4) 
(D.5) 


it  will  be  shown  that 

x(t)  = x0  + xo  (1-exp-t)  + h(t) 

and 

x(t)  = x0exp(-t)  + h(t) 

where  h(t)  and  b(t)  are  random  variables  satisfying 


Eh(t)  - E(h(t)  - 0,  (D.  8) 

E h(t) 2 - %V[2t  - (1-exp-t)  (3-exp- 1)]  (D.9) 

Efi(t)2  - %V[1  - exp-2t]  (D.  10) 

E h(t)h(t)  - %V[1  - exp-t)2  (D. 11) 


(D.  6) 

(D.7)' 


Ev(t)  = 0 

Ev(ta)v(tb)  = V 6(ta  - tb) , 

6(t)  being  the  Dirac  delta  function  satisfying 
6(t)  = 0 if  t 0 and 

oo 

/ 6(t)dt  = 1, 
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and 


« • 


x0  = x(O) , x0  = x(O) 


Proof.  Since 


^ [x(t)expt]  = [x  + x]expt,  equation  (D.l) 


may  be  integrated  to  give 


r t 

x(t)expt  - x0  = / exps  v(s)ds 


Hence 


x(t)  = x0exp(-t)  + h(t)  , 


where 


ft 

h(t)  = J exp(s-t)v(s)ds 


Since 


(D. 12) 


(D. 13) 


(D.  14) 


(D. 15) 


d |*b  . ft 

^ J[  v(s)  [ l-exp(s-t)]  ds  = exp(s-t)v(s)ds , 


it  follows  from  equation  (D. 15)  that 

c ^ 

h(t)  = J v(s)  [l-exp(s-t)]ds 


(D. 16) 


Integration  of  equation  (D. 14)  yields 


x(t)  - x0  + x0  (1-exp-t)  + h(t) 


(D* 17) 


1 
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Equations  (D.6)  and  (D.7)  have  thus  been  derived.  Equations  (D.8) 
easily  follow  from  equations  (D. 15) , (D.16)  and  (D.2).  Equations 
(D. 9)  through  (D. 11)  will  now  be  derived.  From  equations  (D. 15)  , 
(D.16)  and  (D.3)  it  follows  that 


E h(t)2  = V 

•/*  [l- exp  ( s - 1 ) ] 2 ds , 

0 

(D. 18) 

E h(t) 2 = V 

t 

y*[exp2(s- 1)]  ds , 

0 

(D.19) 

E h(t)h(t)  = 

t 

V / exp(s-t)[l-exp(s-t)]ds. 

(D.20) 

Integrating  equation  (D.18)  we  obtain 

E h(t) 2 = V [t  - 2(l-exp-t)  + %(l-exp-2t)] , 

which  after  some  simplification  leads  to  equation  (D.9).  Equations 
(D.19)  and  (D. 10)  are  easily  seen  to  be  equivalent.  From  equation 
(D.20)  we  deduce  that 


E h(t)h(t)  = V[(l-exp-t)  - %(l-exp-2t)] , 


which  can  be  seen  to  reduce  to  equation  (D.ll).  This  completes 
the  derivation  of  the  required  equations. 


Expected  values  of  x,  x,  x2  and  x2  for  large  values  of  t. 

It  follows  from  equations  (D.6)  through  (D.ll)  that 
for  large  values  of  t, 
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Ex(t)  = x0  + Xo 


(D- 21) 


Ex(t)  = 0 


(D. 22) 


Ex(t)2  = (xo  + Xo)2  + Vt  = Vt 


(D. 23) 


Ex(t) 2 = %V 


(D. 24) 


Autocorrelation  functions  for  x and  x. 


If  ta  i t then  we  obtain  similarly  to  equations  (D.18) 


and  (D.19) 


E h(t)h(ta)  = V f [l-exp(s-t)]  [l-exp(s-ta)]  ds  , (D.25) 

o 


E h(t)h(ta)  = V / exp(s-t)exp(s-ta)ds. 


(D. 26) 


Hence, 


E h(t)h(ta)/V  = t-[l-exp-t]-[exp(t-ta)-exp-ta] 


+ %[exp(t- ta)-exp(-t-ta)]  , 


(D. 27) 


E h(t)h(ta)/V  = %[exp(t-ta)-exp(-t- ta)] 


(D. 28) 


Writing  ta  **  t + At, 
we  deduce  that  for  large  t. 


(D. 29) 
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(D. 30) 


E h(t)h(t  + At)  = Vt 

and 

E h(t)h(t  + At)  = %V  exp (-At)  (D. 31)  | 

Thus  for  large  t , 

! 

[Ex(t)x(t  + At)] /[Ex(t)2]=  1 (D.  32) 

and 

[Ex(t)x(t  + At)] / [Ex(t) 2]  = exp(-At)  (D.33) 

A further  quantity  of  interest  is  E[x(t)  - x(t  + At)]  . 

It  follows  from  equations  (D. 6)  and  (D. 8)  that  for  large  t, 

x(t)  - x(t  + At)  = h(t)  - h(t  + At). 

Hence  we  obtain  with  the  aid  of  equations  (D.9),  (D.27)  and  (D.29) 

E[x(t)-x(t  + At)]2  =E.h(t)2  + h(t  + At)2 

- 2 E h(t)h(t  + At) 

= V[t-1 . 5 ] + V[t  + At  - 1.5] 

- 2V[t  - 1 - % exp  (-At)] 

Thus,  for  large  t, 

E[x(t)  -x(t  + At)]2  = V[At  - (1-exp-At)]  (D.34) 

Ji 
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If  furthermore  At  is  small  the  above  may  be  approximated 


E[x(t)  - x(t  + At)]2  = %V  (At)2  (D. 35) 


Change  of  variable  formula  for  the  Dirac  delta  function. 

If  in  equation  (D. 5)  we  change  the  variable  of  integra- 


tion from  t to 


s **  ft 


where  f is  a positive  constant,  we  obtain 


(D. 36) 


<$(s/f)  ds  - f 


(D. 37) 


Consequently, 


6 (s/f ) = f6(s) 


(D. 38) 
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APPENDIX  E 

Selection  of  4 GPS  Satellites  for  Optional  Position  Determination 

In  section  6 of  the  main  text  it  was  shown  that  the 
measurement  cL  to  the  j-th  satellite  is  given  by 

dj  = /(rj-r)T(rj-r)  + b,  (E.l) 

where  r.  is  the  position  vector  of  the  j-th  GPS  satellite,  ? is 
the  position  vector  of  the  user,  and  b is  a measurement  bias. 

In  this  appendix  we  shall  determine  where  the  4 GPS 
satellites  ideally  should  be  located  in  order  that  the  user's 
position  may  be  calculated  with  the  least  amount  of  error.  To 
do  that  we  make  the  following  assumptions 

E fir.  = E fid.  = 0,  (E.  2) 

J J 

where  fir^  is  the  error  in  the  position  of  the  j-th  GPS  satellite 
and  fidj  is  the  measurement  error.  We  further  assume  that  the 
errors  are  uncorrelated: 

E fir^fir^  = 0,  E fid^fid^  = 0 for  j ^ k, 

and 

E fir^fir..1  = o2  I,  E(5d..)2  = c02  (E.3) 

where  o and  Oo  are  scalars,  and  I is  the  3x3  identity  matrix. 
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If  the  error  in  the  calculated  user  position  is  denoted 
by  6r  and  the  error  in  the  calculated  bias  by  5b,  then  it  follows 


from  equation  (E.l)  that  if  the  errors  are  small  then 


6d.  = v.1  (6r.  - 6r)  + 6b, 

J j J 


where  the  unit  vector  v^  is  defined  by 


vj  ■ 


Rearranging  the  terms  in  equation  (E.4)  we  find  that 


v.  6r  - 6b  =»  s . , 
J J 


where 


8j  - 


(E.4) 


(E.5) 


(E.6) 


(E.7) 


From  equations  (E.2)  and  (E.3)  we  deduce  that 


E Sj  = 0, 


(E.8) 


Esisk  = 0 if  j y4  k 


(E'.  9) 


E Sj 2 = v/  o2  Vj  + a o: 


Since  v^  is  a unit  vector  the  last  equation  reduces  to 


E Sj 2 - o2. 


(E.  10) 
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where 


o2  = a2  + oo2 


(E. 11) 


Combining  the  measurement  from  four  GPS  satellites,  equation 
(E.6)  may  be  rewritten  as 


[«] 


= s , 


where 


r t ,i 

V)  , -1 

and  s = 

sf 

T 

V2  . -1 

s2 

V?  . -1 

s3 

T 

vs,  , -1 
■»  « 

S4 

(E.  12) 


(E. 13) 


It  is  interesting  to  note  from  equations  (E.12)  and 
(E.13)  that  the  position  error  6r  is  a function  of  the  'user  to 
GPS  satellite'  direction,  but  independent  of  the  corresponding 
distance. 

From  equations  (E.8),  (E.9),  (E.10),  and  (E.13)  we 
deduce  that 


Es  3 0 and  Ess  = a2  I, 

where  I is  the  4x4  identity  matrix. 

Solving  equation  (E.12)  we  obtain 


(E.  14) 


□ 


= A”  1 s 


(E.  15) 
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We  hence  deduce  with  the  aid  of  equation  (E.14)  that 
E 

and 


It  will  now  be  shown  that  in  order  that  the  expected 

_ _ *T» 

square  of  the  position  error  (trace  E 6r6r  ) be  a minimum,  the 
angles  between  the  four  vectors  vi,  V2,  v3  and  v4  should  all  be 
equal. 

Let 


> 

1 

M 

II 

r»i 

a= 

i 

H 

X 

1 

bj  b2  b3  b4 
hi  h2  h3  h4 


(E.  18) 


where  B is  a 3 x 4 matrix;  h is  a 4- vector;  bi,  b2,  b3  and  b4  are 
3-vectors;  hi,  h2,  h3  and  h4  are  scalars. 

Since  AA" 1 = A~lA  = I it  follows  from  equation  (E.13) 
and  (E. 18)  that 


v±  bj  - Slj  + hj 


where  6^j  is  the  Kronecker  delta. 


(E. 19) 


2 b v.J 
i=l  1 1 


I, 


(E.  20) 


1 


6r 

= 0 

(E. 16)  j 

6b 

! ; 

|*6r  6rT  6r6b 

= o2  A-1A-T  (E.  17) 

| 

6b 

5rT  6b2 

L 
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2 b.  = 0, 
i=l  1 


(E. 21) 


2 h.  = -1 
i=l  1 


(E. 22) 


From  Equations  (E.17)  and  (E.18)  we  deduce  that 


the  Lagrangian  multipliers  A.^  having  been  introduced  to  take  into 
account  the  fact  that  each  vi  is  a unit  vector.  The  following 
result  is  derived  at  the  end  of  this  appendix 


- trace  BBT  = -2  BBT  b. 
3vj  J 


(E. 25) 


Using  equation  (E.25)  we  deduce  from  equation  (E.24) 


that 


A.v.  = BBt  b. 

J J 


(E. 26) 


It  follows  from  equations  (E.21)  and  (E.26)  that 


2 A.v.  = 0 


(E. 27) 
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trace  E 6r6rT  = cr2  trace  BBT 


(E. 23) 


In  order  that  the  above  quantity  be  a minimum  it  is 


necessary  that 


•3 — trace  BB1  + 2 (v.  v.  - 1)A.  = 0 , 
3v.  i=l  1 1 x 


(E.24) 


T 

Premultiplying  equation  (E.26)  by  we  deduce  with  the  aid 
of  equation  (E.19)  that 


BBT  bj  = A^  bj 


Vj 


" Xj(6ij 


+ hi> 


(E.  28) 


Since  the  left  hand  side  of  the  above  equation  is  symmetric  in 
i and  j we  conclude  that 


Xj  hi  “ Xi  hj 


(E. 29) 


Summing  the  above  equation  with  respect  to  i we  deduce  with  the 
aid  of  equation  (E.22)  that 


4 

2 

i=l 


Xi 


(E.  30) 


There  does  not  appear  to  be  any  reason  why  A^  should  differ  from 
A^.  We  therefore  make  the  assumption  (later  to  be  justified,  of 
course)  that  for  all  j , 


A 


j 


= A 


It  hence  follows  from  equation  (E.30)  that 


Since  by  equation 


hj  * 

(E.  18) 


BB  - Z b.b,T 
j = l J J 


(E.  31) 


(E.  32) 


(E.  33) 
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we  conclude  from  equations  (. 31) , (E.20)  and  (E.26)  that 


XI  = BBT  BBT  . 


(E. 34) 


Since  BB1  is  a semi-positive  definitive  matrix  this  is  only 
possible  if 


bbt  = /X  I 


(E. 35) 


It  hence  follows  from  equation  (E.26)  that 


bj  = /X  Vj  . 


(E. 36) 


From  the  above  and  equation  (E.32)  and  (E.19)  we  conclude  that 


/X  Vi  v3  - 4jJ  - 4 


(E.  37) 


Since  v^^  is  a unit  vector  it  follows  that  (i  = j in  the  above 
equation) 


/A  « ■$  . 


(E.  38) 


Hence  if  i f j , 


T j. 

vi  vj  - - * 


(E.  39) 


This  is  the  desired  result.  We  conclude  from  equations 
(E. 17) , (E.18),  (E. 35) , (E.38),  (E.32)  and  (E.21)  that 


E Ur6xT  6r<5b 


<5b<5rT  6b2 


(E. 40) 
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A set  of  unit  vectors  satisfying  equation  (E.39)  are 

given  by 


T 

Vi 

= ( 

0 

, 0 

1 

) 

T 

Va1 

- ( 

2/2/3 

, 0 

, -1/3 

) 

T 

V31 

= ( 

-/in 

. /Tn 

, -1/3 

) 

T 

V41 

= ( 

-/in 

, -/in 

, -1/3 

) 

By 

equations 

(E. 36)  and 

(E.38) 

(E.41) 


(E. 42) 


With  hj  given  by  equation  (E.32),  equations  (E.19) 
through  (E.22)  are  readily  verified. 

From  equations  (E.33),  (E.20)  and  (E.42)  it  follows  that 

BBt  - $ I,  (E, 43) 


which  is  consistent  with  equations  (E. 26) , (E.42),  (E. 31)  and 
(E.38).  The  solution  set  (E.41)  is  thus  justified.  Note  that  the 
angles  between  the  vectors  v^  equal  cos-1  (-1/3)  or  109.5  degrees, 
and  that  the  expected  square  of  the  position  error  as  given  by 
equations  (E.23)  and  (E.43) 

E <5rT6r  = * o2  (E.44) 
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E.  1 


The  Selection  of  Two  Satellites  when  two  have  already 
been  Selected. 

Let  us  assume  that  we  are  given  two  vectors  Vi  and  V2 . 
There  is  no  loss  of  generality  in  assuming  that  they  are  of  the 
form 


VjT=  (cosp,  sinp,  o) 
v2t=  (cosp, -sinp,  o) 


(E.  45) 


For  reasons  of  symmetry  it  follows  that  v3  and  v4  must  be  of  the 
form 


v3t  = (cosq,  0 , sinq) 
v,Ts=  (cosq,  0 , -sinq), 


(E.  46) 


for  some  angle  q to  be  determined.  From  the  above  and  equation 
(E.13)  it  follows  that 


cosp , 

sinp, 

0 , 

-1 

cosp. 

-sinp. 

0 , 

-1 

cosq. 

0 . 

sinq  , 

-1 

cosq, 

0 

-sinq  , 

-1 

(E.  47) 


1 


Inverting  the  above  matrix  we  obtain 
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T 

Since  BB  is  proportional  to  the  expected  square  of  the 
position  error,  we  choose  q such  as  to  minimize  that  trace.  Letting 
the  partial  derivative  with  respect  to  q vanish  we  deduce  that 


2sing + = 0 (E.51) 

(cosp-cosq)3  sin3q 

Let  cosp  = c and  cosq  = x (E.52) 

Corresponding  to  equation  (E.51)  we  then  obtain 

f(x,c)  - 0,  (E. 53) 

where 

f(x,c)  - 2(l-x2) 2 + x (c-x) 3 (E. 54) 

In  finding  the  solutions  of  the  above  equations,  we  may  assume  that 

c > 0,  (E. 55) 


since  this  only  involves  the  definition  of  the  coordinate  axes. 
Obviously,  c < 1.  From  equation  (E.54)  we  find  that 

f (~°°,  c)  = + °° 
f(-l.c)  » - (c+1)  3 
f(0,c)  = 2 
f(l,c)  = -(1-c)3 

f (°°,  c)  = + CO 
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From  the  above  it  is  evident  that  equation  (E.53) 
always  has  4 real  solutions , and  of  those  there  is  always  one 
and  only  one  in  the  interval  (-1,0).  Since  cosp  is  positive  it 
is  quite  obvious  from  equation  (E. 50)  that  the  desired  solution 
is  negative.  Although  it  is  not  simple  to  obtain  an  exact 
solution,  a good  approximation  is  given  by 


where 


x = g(c). 


g (c)  = .327c  - .765 


(E. 56) 


(E. 57) 


the  linear  approximation  being  based  on  the  exact  solutions  for 
c = 1 [x  = (1/17- 5)/ 2]  and  c = 0 [x  = V2-/2]. 

Let  A2  denote  the  expected  square  of  the  position  error. 


From  equations  (E.23),  (E. 50)  and  (E.52)  we  obtain. 


A 2/a2  = 


(c-x)2  2 (1-c) 2 2 (1-x2) 


(E.  58) 


A comparison  of  exact  and  approximate  values  of  x as 
functions  of  c are  given  in  table  E.l.  Formula  (E.58)  has  also  been 
evaluated  in  the  table.  Note  that  the  optimum  value  for 
c =*.577  = 1//)  corresponds  to  solution  (E.41). 
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TABLE  E.l 


Exact  and  Approximate  Values  of  x,  and  the 


Expected  Square  of  the  Position  Error  as  a 
Function  of  c. 


c | x(exact)  s 

| x = g(c) 

9 1 2 

0 3 4 

42/o2 

.000 

-.765 

-.765 

180.0 

80.2 

3.41 

.174 

-.710 

-.708 

160.0 

89.9 

2.80 

.342 

-.655 

r.653 

140.0 

98.5 

2.45 

.500 

-.603 

-.602 

120.0 

106.0 

2.27 

.643 

-.555 

-.555 

100.0 

112.6 

2.27 

.766 

-.515 

-.515 

80.0 

118.0 

2.50 

.866 

-.482 

-.482 

60.0 

122.4 

3.20 

.940 

-.458 

-.458 

40.0 

125.5 

5.44 

.985 

-.443 

-.443 

20.0 

127.4 

17.91 

1.000 

-.438 

-.438 

0.0 

128.0 

CO 

,577 

-.577 

-.576 

109.5 

109.5 

2.25 

0i2  is  the  angle  between  satellites  1 and  2 
as  seen  from  the  user. 

03i,  is  the  angle  between  satellites  3 and  4 
as  seen  from  the  user. 

A2  is  the  expected  error  of  the  square  of 
the  position  error 

o2  is  the  sum  of  the  measurement  error  variance 
and  the  variance  of  the  GPS  position  error 
measured  along  any  axis. 
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(E.  60) 


(E. 61) 


(E. 62) 
(E.  63) 
(E.  64) 
(E. 65) 


(E.  66) 


Fro©  the  above  and  equation  (E.65)  we  find  that 


hi,  = (1  - vj  Ao 'dj)"1 


(E.  67) 


From  equation  (E.62)  and  (E.64)  we  obtain 


Hence , 


v*  (Ao 1 + Ao 1 doho ) - ho  = 0. 


h?  = (1-vi?  A^Mo)"1  vf  Ao1 


= hi*Vi,T  Ao1 


(E. 68) 


From  the  above  and  equation  (E.62)  we  deduce  that 


B0  = AJ*  + AoldoV?  Ao*hi, 


(E. 69) 


Writing 


e o = Af 1 d0  and  f = AJ x Vi,  , 


(E.  70) 


equations  (E.66),  (E.67)  and  (E.69)  may  be  rewritten  in  the  form 


and 


b 4 = - e ohi, 

hi,  = (l-vfe0)_1 
B0  = Ao"1  + e0  fT  hi, 


(E. 71) 
(E. 72) 
(E. 73) 


Comparing  equation  (E.18)  and  (E.61)  we  see  that 


B - [Bo  , bH] 


(E.  74) 
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It  hence  follows  from  equations  (E.71)  and  (E.73)  that 


BBt  = Ao 1 Ao  ^ + Ao 1 f e0T  hi,  + e 0 f Ajj^  hi, 


+ e0fTf  ej  hi;  + e0e?  h* 


(E. 75) 


Hence 


trace  BBT  = trace  Ao^Aj;7  + 2hi,e?Ao1f 


+ h*  e0  e0  (1  + f f) 


(E. 76) 


Since  the  expected  square  of  the  position  error  is  proportional 
T 

to  trace  BB  , it  follows  that  we  must  choose  the  fourth  satellite 
such  that 

2h„e0T  Ao1  f + K e^e0  (1  + fT  f) 

is  a minimum.  Ao  as  defined  by  equation  (E.59)  may  be  inverted 
using  the  formula 


Ao 1 = (v2Av3,  v3AVi,  ViAv2)/ [(viAv2) -v3]  (E.77) 

E. 3 Derivation  of  Equation  (E.25) 

The  desired  formula  is  most  easily  derived  using  the 
customary  index  summation  convention.  In  what  follows  Latin  indices 
will  assumethe  values  1,2, 3, 4 and  Greek  indices  the  values  1,2,3. 
Since  B is  a 3 x 4 matrix  it  follows  that 


<BB  >aB  - BaiBBi 


(E.  78) 


Hence 


trace  (BBT)  = BaiBai 


(E. 79) 


-85- 


(E.  80) 


From  equation  (E. 18)  we  deduce  that 


<*W«B  ■ SaB 


Consequently 


and 


BakAks  ^ots  * 


A.  <B»kAks>  - 0 
jy 


(E. 81) 


(E. 82) 


Hence 


6 =0 
sy 


(ax.  B«k)Aks  + Vj1 

J r* 

Post  multiplying  the  above  equation  by  Asi  we  deduce  that 


(E. 83) 


1_  B . + B .A"\  = 0 
8Aj^  ax  aj  yr 


x.  e. 


-B  .B  . 
aj  yi 


Bod  iJL,  Bai  “ BajByiB«i' 


and 


JK 


JU 


3 trace  BBT  = -2  B^B^B^ 


jy 


(E . 84) 

(E. 85) 

(E. 86) 
(E. 87) 


In  accordance  with  equations  (E.13)  and  (E.18)  this  may  also  be 
written  as 

(k  trace  BbT)m  - • 2 <BBT)M«<Va  (E'88) 
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Hence 


trace  BB^  = -2BB  b ^ , 


the  desired  result. 
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(E.25) 


appendix  f 

Flow  Charts  and  Short  Descriptions  of  Fotonap  Subroutines 


>C1HZH 


o w 


pdWHHJOfi 


ALPHABETIC  LISTING  AND  SHORT  DESCRIPTION  OF 


PHOTONAP  SUBROUTINES 

SUBROUTINE 


Name 

CALLED 

BY 

PURPOSE  OR  FUNCTION 

ATMIN 

FINDOG 

Reads  tabulated  atmospheric  densities 

from  file  LUA  1 

ATMOUT 

GENFIL 

! 

Writes  tabulated  atmospheric  densities 
on  file  LUA 

BETA 

ION 

Used  in  computation  of  ionospheric 
corrections 

CAT000 

SPOLCD 

Converts  "free  form"  card  input  to 
standard  NAP  format  (not  on  1108) 

DARITE 

SOLVER 

INVPRE 

REDS KM 
INVERT 
INVRTB 
INVRTC 

Utility  routine  for  direct  access 
file  43 

DAYHMS 

NUTION 

EDIT 

INTPRT 

IONOSF 

RESID 

SIMOUT 

VISIBI 

KMNOUT 

KMNSIM 

KMNSM2 

Converts  (Julian  day,  seconds  of  day) 
to  (year,  month,  day,  hour,  min,  sec) 

DBREAD 

PREINT 

Utility  routine  for  direct  access 
file  41 

DBURN 

INTEG 

. 

Adds  velocity  increment  to  satellite 
velocity  (Discrete  thrust) 

DEFALT 

DREDIT 

Initializes  program  constants  to  their 
default  values  (See  User's  Guide) 

DENS 

EXPAND 

Computes  atmospheric  density  as 
function  of  height  (*) 

DGRITE 

PTSTAR 

SOLVER 

PRTIAL 

Utility  routine  for  direct  access 
file  40 
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A 


SUBROUTINE 


1 


Name 

CALLED 

BY 

PURPOSE  OR  FUNCTION 

DKGK 

ION 

Used  in  computation  of  ionospheric 
corrections 

DKSICO 

ION 

Used  in  computation  of  ionospheric 
corrections 

DRAGU 

INTEG 

Called  at  the  end  of  each  drag 
segment  (except  the  last)  to  calcu- 
late the  appropriate  partial  deriva- 
tives 

DRAVAR 

ENGRAT 

Computes  drag  contribution  to  varia- 
tional equations  (*) 

DREDIT 

OLDMAN 

Control  routine  for  processing  control 
card  input 

DUMDUM 

PHOTO 

Dummy  routine  used  for  switching 
program  overlays 

EDIT 

DREDIT 

Edits  observed  data  or  (in  simulation 
mode)  generates  random  numbers 

EIGN 

FINALP 

Computes  eigen-values  of  a matrix 

ENBVAR 

ENGRAT 

Computes  central  term  and  planetary 
contribution  to  variational  equations 
(*) 

ENGRAT 

INTEG 

Control  routine  for  each  integration 
step  (*) 

ENROOT 

ENGRAT 

OCCULT 

Finds  the  zero  of  a function  expressed 
as  a power  series 

EXPAND 

ENGRAT 

Develops  power  series  coefficients  for 
satellite  vector  (*) 

FINALP 

OLDMAN 

Prints  final  results 

FINDOG 

INTEG 

Initializes  integrator  at  start  of 
integration  or  on  change  of  origin 

FIREAD 

FINDOG 

INTGA 

Utility  routine  for  direct  access 
file  41 

FLREAD 

FINALP 

Utility  routine  for  direct  access 
file  41 

J 
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. 

. 

SUBROUTINE 

CALLED 

Name 

BY 

PURPOSE  OR  FUNCTION 

Debug  print  for  subroutine  GENFIL 
(formerly  part  of  GENFIL) 

Generates  internal  files  based  on 
control  card  input 

Routine  for  handling  geoceiver 
measurements 

Reads  data  tape  in  GEOS  format 

Reads  planetary  ephemerides  from 
file  (10) 

Used  in  computation  of  ionospheric 
corrections 

Reads  data  tape  in  NAP  format 

Estimates  ground  point  coordinates  by 
projecting  photographic  plate  coordi- 
nates onto  Earth’s  surface 

Tropospheric  refraction  corrections 
for  geoceiver  data 

Subroutine  for  unpacking  integers 
and  storing  them  (unpacked)  in  an 
array 

Function  for  packing  integers 

Function  of  unpacking  integers 


Called  by  INPCRD  for  debug  print 

Processes  control  card  input 

Used  for  processing  series  100  input 
cards  (formerly  part  of  routine 
INPCRD) 
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SUBROUTINE 


Name 


INP200 


INP300 


INP600 


INP700 


INTCOD 


INTEG 

INTERP 


INTPRT 


INVERT 


IN VP RE 


INVRTB 


CALLED 

BY 


INPCRD 


INPCRD 


INPCRD 


INPCRD 


INP200 


OLDMAN 

TIMARR 


ENGRAT 

INTGA 

PREPAR 

ENGRAT 

PRTIAL 

PTSTAR 

NUTION 


SOLVER 

SOLVER 


SOLVER 


PURPOSE  OR  FUNCTION 


Used  for  processing  series  200  input 
cards  (formerly  part  of  routine  INPCRD) 

Used  for  processing  series  300  input 
cards  (formerly  part  of  subroutine 
INPCRD) 

Used  for  processing  series  600  input 
cards  (formerly  part  of  routine  INPCRD) 

Used  for  processing  series  700  input 
cards  (formerly  part  of  routine  INPCRD) 

Generates  arrays  for  recovery  of  gravity 
parameters  (spherical  harmonics  or 
mascons) 

Control  routine  for  integrator 

Interpolation  routine  used  for  setting 
up  tables  of  time  corrections 

Prints  integrator  output  at  end  of 
integration  of  each  arc 

Prints  time  corresponding  to  integrator 
output 

Estimates  difference  between  integrator 
time  and  UTC  through  interpolation 


Estimates  difference  between  UT1  and 
integrator  time  through  interpolation 

Part  1 of  matrix  inversion 

Used  in  processing  photogrammetric  data, 
(i)  resequences  file  26  (IPASOT)  of 
ground  point  images  on  file  38.  (ii) 
on  first  iteration  writes  ground  point 
records  from  file  25  (IARCOT)  to  file 
23  (IGPT) 

Part  1 of  matrix  inversion.  Prints 
intermediate  results 
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SUBROUTINE 


Name 

CALLED 

BY 

PURPOSE  OR  FUNCTION 

INVRTC 

SOLVER 

Similar  to  INVRTB,  but  solution 
obtained  without  completing  matrix 
inversion 

INVSYM 

KMNCON 

IN  SITU  inversion  routine  for  positive 
definite  matrix  in  upper  diagonal  form 

INVSYS 

EDIT 

Similar  to  INVSYM.  Additional  feature 
is  to  check  for  singularity 

ION 

IONOSF 

Used  in  computation  of  ionospheric 
corrections 

IONOSF 

PRTIAL 

Used  in  computation  of  ionospheric 
corrections 

JACHIA 

DENS 

Computes  atmospheric  density  (Lockheed- 
Jacchia  model) 

JULDAY 

SPOLCD 

INP200 

IONOSF 

Converts  (year,  month,  day,  hours, 
minutes,  seconds)  to  (Julian  day, 
seconds  of  day) 

KEPLER 

PREINT 

Converts  Keplerian  to  Cartesian  input 

KICKER 

INTEG 

Initializes  integrator  common  blocks 

KMNCON 

OLDMAN 

Control  routine  for  Kalman  filtering 
and  smoothing 

KMNEVA 

KMNCON 

Used  in  Kalman  filtering  for  evaluating 
integrated  power  series 

KMNIDE 

KMNINI 

KMNCON 

Store  (9x9)  identity  matrix  in  required 
location 

KMNINI 

KMNCON 

Initialization  routine  for  Kalman  filtering 
and  smoothing 

KMNINV 

KMNCON 

Used  in  Kalman  filtering  for  calculating 
the  transition  matrix  relative  to  the 
previous  time  point  (given  the  transition 
matrices  relative  to  epoch) 

KMNMP1 

KMNCON 

T 

Utility  routine  for  computing  C = A * B 
where  A is  symmetric  and  stored  in  upper 
triangular  form 
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SUBROUTINE 


Name 


KMNMP2 


KMNMP3 


KMNOUT 


KMNRAN 


KMNSIM 


KMNSM2 


MAGFIN 


MAIN 

(see  photo) 
MASCON 


MATINV 


MATZEV 


MATZEX 


CALLED 

BY 


KMNCON 


KMNCON 


KMNCON 


KMNCON 


KMNCON 


KMNCON 


EXPAND 


KMNINV 

KMNCON 

DBURN 

DRAGU 

ENGRAT 


ENGRAT 


PURPOSE  OR  FUNCTION 


Utility  routine  for  computing  C = A * B 
where  C is  symmetric.  C is  stored  in 
upper  triangular  form 

Utility  routine  for  computing 
C = AT  * b * A,  and  x = A^y.  B and  C 
are  symmetric  and  stored  in  upper 
triangular  form 

Used  to  output  the  parameter  estimates 
(and  their  covariances)  from  the  Kalman 
filter  and  smoother.  (Output  on  file 
29  (ITAPE)) 

Random  number  generator  associated 
with  simulations  in  Kalman  filtering 

Simulates  GPS  output  used  in  Kalman 
filtering  (Output  on  file  33  (LF33)) 

Printout  routine  for  Kalman  filtering. 
In  simulation  mode  prints  full  9 para- 
meter state-vector.  In  filter  or 
smoother  mode  prints  only  3 parameter 
state-vector  (first  6 parameters 
printed  in  KMNOUT) 

Used  in  computation  of  ionospheric 
corrections 

Control  routine  for  NAP  program 


Computes  mascon  contributions  to 
satellite  acceleration  (**) 

Matrix  inversion 


Evaluates  state  transition  matrix  at 
end  of  integration  step  (*) 

Develops  coefficients  for  power  series 
expansion  of  state  transition  matrix 


SUBROUTINE 


Name 


MAT3T3 


MEASUR 


MEASXX 


NBDNEX 


NBDPEX 


CALLED 

BY 


MEASUR 

MEASXX 

PLATEC 

GROUND 

PTSTAR 

TERRA 

MEASUR 

SP02A3 

MEASXX 

TDRSS 

PTSTAR 

PLATEC 

PTSTAR 

NUTION 

PTSTAR 

PLATEC 

STARAN 

ROTATE 

KMNCON 


PRTIAL 


PRTIAL 


ENGRAT 


ENGRAT 

EXPAND 


PURPOSE  OR  FUNCTION 


PLATEC  I Matrix  multiplication  (2x3)x(3x3) 


Matrix  multiplication  (3x3)'x(3x3) 
('  indicates  transpose) 


Matrix  multiplication  (3x3)x(3xl) 


Matrix  multiplication  (3x3)x(3x2) 
Matrix  multiplication  (3x3)x(3x3) 


C = A * B.  A is  a (9x9)  matrix, 

C and  B are  (9xn) , where  n = 9 (entry 
MAT99)  or  n = 1 (entry  MAT91) 

Routine  for  handling  the  following 
measurement  types:  RANGE,  AZIMUTH, 
ELEVATION,  RIGHT  ASCENSION,  DECLINATION, 
MINITRACK  DIRECTION  COSINES , X30  and  Y30 
ANGLES,  DISTANCE  TO  ELLIPSOID,  RANGE 
RATE,  MINITRACK  RATES,  X85  and  Y85 
ANGLES,  STATE  VECTOR  MEASUREMENTS 

Routine  for  handling  the  following 
measurement  types:  RANGE  SUM,  RANGE 

SUM  RATE,  GRARR,  TDRSS,  GEOCEIVER 

Develops  coefficients  for  power  series 
expansion  of  Sun,  Moon,  and  Planets  (*) 

Computes  central  term  and  planetary 
contribution  to  satellite  acceleration 


SUBROUTINE 


1 


I 


SUBROUTIN 

Name 

E 

CALLED 

BY 

' 

PURPOSE  OR  FUNCTION 

NEWJPL 

INTEG 

Used  in  conjunction  with  READE  and 

GETTAP  to  obtain  planetary  ephemerides 

NUTION 

GENFIL 

Calculates  precession/nutation  matrix 

PREINT 

and  Greenwich  Hour  Angle  j 

I 

, 

ENGRAT 

PRTIAL 

PTSTAR 

SOFRT 

KMNCON 

OCCULT 

ENGRAT 

If  satellite  is  orbit  around  body  A, 

FINDOG 

this  routine  determines  if  satellite 
is  visible  from  body  B 

OFDATE 

SOFORT 

Rotates  vector  or  matrix  from  "inertial 

SP02A3 

1950.0"  to  "true  of  date"  (Double 

SOFSEC 

XFORM 

KMNEVA 

Precision) 

OLDMAN 

PHOTO 

("old  main")  secondary  control  routine 
for  NAP  program 

04DATE 

SOFRT 

Rotates  vector  or  matrix  from  "inertial 

1950.0"  to  "true  of  date"  (single 
precision) 

PAGE 

SPOLCD 

DREDIT 

EDIT 

INTCOD 

PREPAR 

RESID 

FINALP 

SIMOUT 

VISIBI 

KMNINI 

SOLVER 

SOLV2 

INVRTB 

INVRTC 

Prints  page  heading 

PFSOLV 

SOFORT 

Evaluates  partials  of  satellite  vector 
w.r.t.  continuous  thrust  parameters 
using  previously  computed  power  series 
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SUBROUTINE 


Name 

CALLED 

BY 

PURPOSE  OR  FUNCTION 

PFVARY 

ENG RAT 

j 

Develops  power  series  coefficients  for 
partials  of  satellite  vector  w.r.t. 
continuous  thrust  parameters  (*) 

PLATEC 

PTSTAR 

Computes  predicted  photographic  plate 
coordinates  and  range  (Photogrammetric 
measurement  types  7-9) 

POTSOL 

SOFORT 

SOFSEC 

Evaluates  partials  of  satellite  vector 
w.r.t.  gravity  parameters  using  pre- 
viously computed  power  series 

POTVAR 

ENGRAT 

Computes  central  body  (excluding  central 
term- -see  ENBVAR)  contribution  to 
variational  equations  (*) 

PREINT 

PREPAR 

Sets  up  arrays  for  integrator  based 
on  current  values  of  "solve  for" 
parameters 

PREPAR 

OLDMAN 

Sets  up  arrays  for  integrator  based 
on  control  files 

PRTBG1 

PRTIAL 

Output  debug  print  from  subroutine 

PRTIAL 

PRTBG2 

PRTIAL 

Output  debug  print  from  subroutine 

PRTIAL 

PRTBG3 

PRTIAL 

Output  debug  print  from  subroutine 

PRTIAL  1 

PRTIAL 

OLDMAN 

Computes  differences  between  observa- 
tions and  predicted  observations.  Also 
computes  associated  partials.  Results 
output  on  file  (ISFILE) 

PTCMPA 

PRTIAL 

Data  compression  associated  with  sub- 
routine PRTIAL  (formerly  part  of  PRTIAL) 

PTINIT 

PRTIAL 

Used  for  initializing  variables  used  in 
subroutine  PRTIAL  (formerly  part  of 

PRTIAL) 

PTSTAR 

PRTIAL 

Routine  for  handling  photogrammetric 
measurements  (formerly  part  of  PARTIAL) 
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f 

t 
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SUBROUTINE 


Name 

CALLED 

BY 

! 

PURPOSE  OR  FUNCTION 

RAN601 

INP600 

INTCOD 

1 

Modifies  a parameter  value  by  adding 
a random  Gaussian  number  with  a given 
standard  deviation 

READE 

NEW  J PL 

Used  in  conjunction  with  GETTAP  to 
obtain  planetary  ephemerides 

READER 

EDIT 

RECPOT 

PRTIAL 

SOFORT 

SOFSEC 

PTSTAR 

Utility  routine  for  read/write  from/to 
sequential  file 

READU 

ION 

Used  in  computation  of  ionospheric 
corrections 

RECCOF 

ENGRAT 

Develops  power  series  coefficients  for 
partials  of  satellite  w.r.t.  a single 
parameter  (Used  for  solar  pressure  and 
drag)  (*) 

RECPOT 

ENGRAT 

Develops  power  series  coefficients  for 
partials  of  satellite  vector  w.r.t. 
gravity  parameters  (*) 

REDISK 

PRTIAL 

PTINIT 

PTSTAR 

j 

Utility  routine  for  read/write  of  \ 
totally  stable  parameters  on  random  1 
access  file  \ 

REDSKK 

SOLVER 

SOLV1 

INVERT 

INVRTB 

INVRTC 

Utility  routine  for  read/write  of  totally  j 

stable  parameters  on  random  access  file  j 

REDSKM 

SOLVER 

SOLVFU 

SOL  VI 

SOLV2 

INVRTB 

Utility  routine  for  read/write  of  normal 
equation  coefficient  matrix  on  random 
access  file 

REFRCT 

PRTIAL 

TDRSS 

Computes  tropospheric  refraction  corrections 

REPRT2 

DREDIT 

Generates  printed  report  of  run  conditions 
as  specified  by  the  control  cards 

(temporarily  removed  from  NAP) 

-109- 

. . IA 


SUBROUTINE 


Name 

REREAD 


CALLED 

BY 


RESID 


PURPOSE  OR  FUNCTION 

Utility  routine  for  direct  access 


file  41 


RESID 


OLDMAN 

FINALP 


Prints  measurement  residuals 


ROTATE 


PTSTAR 

STARAN 

TERRA 


Computes  a rotation  matrix  corresponding 
to  sequential  rotations  about  principal 
axes 


ROTINT 


ROTPAR 


PRTIAL 


MEASUR 

MEASXX 


Converts  (latitude,  longitude,  height) 
to  Cartesian  coordinates.  Computes 
rotation  matrix  "Earth  fixed  Geocentric 
to  Local" 

Rotates  measurement  partials  w.r.t. 
satellite  state-vector  from  "Earth  fixed" 
to  "True  of  date" 


ROTVFD 

ROT1 

RSUM 

SELECT 

SICOJT 

SIGWT 


MEASUR 

MEASXX 

SOFSEC 

PTSTAR 

ROTATE 

MEASXX 


KMNCON 


ION 

MAGFIN 

FINALP 


Rotates  a vector  (v)  from  "true  of  date" 
(D)  to  "Earth  fixed"  (F) 


Computes  a rotation  matrix  corresponding 
to  a rotation  about  a principal  axis 

Computes  predicted  range  sum  and  range 
sum  rate  measurements  (Measurement 
types  16-17) 

Associated  with  GPS  measurements.  Computes 
GPS  satellite  position,  user  distance 
to  them  and  partial  derivatives  w.r.t. 
user  position.  In  simulation  mode, 
selects  an  optimal  set  of  4 GPS  satellites. 

Used  in  computation  of  ionospheric 
corrections 

Converts  weights  to  standard  deviations 
and  vice  versa 


SIMOUT 


OLDMAN 


Computes  simulated  measurements.  Outputs 
results  on  file 


SKRIV 


PRTIAL  Utility  routine  for  writing  data  on 

PTSTAR  sequential  file 


-110- 


SUBROUTINE 


Name 

CALLED 

BY 

PURPOSE  OR  FUNCTION 

SOFORT 

PRTIAL 

MEASUR 

GEOCEG 

PTSTAR 

Reads  power  series  for  primary  satellite 
state  vector  and  partials  from  sequential 
file  and  evaluates  at  required  time 
point 

SOFRT 

FIN ALP 

Reads  power  series  for  primary  satellite 
state-vector  and  partials  from  sequential 
file  and  evaluates  partials  w.r.t. 
the  initial  state-vector . Computes 
state-vector  covariance  matrix  1 

SOFSEC 

MEASXX 

GEOCEG 

Reads  power  series  for  secondary  satellite 
state-vector  and  partials  from  sequential 
file  and  evaluates  at  required  time  point 

SOLAWA 

SOLVFU 

Function  for  computing  A^WA,  ATWy,  yTWy. 

A is  a (6xn)  matrix,  y is  a 6-vector. 

W is  a (6x6)  symmetric  matrix  stored  in 
upper  triangular  form.  (called  from 
subroutine  SOLVER  when  processing  data 
from  Kalman  filter  output) 

SOLREC 

SOFORT 

SOFSEC 

Evaluates  power  series  to  obtain 
partials  of  satellite  state-vector 
w.r.t.  a parameter  (solar  pressure  and. 
drag) 

SOLVER 

OLDMAN 

Control  routine  for  generating  and 
solving  normal  equations 

SOLVFU 

SOLVER 

Used  for  computing  the  contribution  of 

6 correlated  measurements  to  the  Normal 

Equations  Matrix  and  Vector 

SOLV1 

SOLVER 

Used  for  computing  contribution  of  a 
priori  parameter  values  to  Normal 

Equations  Matrix  and  Vector  (Subroutine 

S0LV1  was  formerly  part  of  subroutine 

SOLVER) 

SOLV2 

SOLVER 

Prints  correlation  vector  for  each  arc 
and  stores  primary  arc  covariance  matrix 
(entry  S0LV2) . Stores  parameter  numbers 
and  stability  types  for  primary  and 
secondary  arcs  (entry  S0LV2A)  Initializes 
normal  equations  matrix  and  vector  for 
multiple  drag  segments  (entry  S0LV2B). 

(Subroutine  S0LV2  was  formerly  part  of 

subroutine  SOLVER) 
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SUBROUTINE 


Name 

SOLV3 

SOREAD 

SORITE 

SPOLCD 

SPOOL 1 

SPOOL2 

SP0100 

SP02A3 


CALLED 

BY PURPOSE  OR  FUNCTION 


SOLVER 


Used  for  deciding  when  to  terminate 
iterative  solution  process  (entry 
S0LV3) . Prints  summary  associated 
with  inversion  of  Normal  Equations 
Matrix  (entry  S0LV3A)  (Subroutine 
S0LV3  was  formerly  part  of  subroutine 
SOLVER) 


SOLVER 


Utility  routine  for  reading  a sequential 
file 


SOLVER 


Utility  routine  for  writing  a sequential 
file 


OLDMAN 


SPOLCD 


SPOLCD 


SPOLCD 


MEASUR 


Scans  NAP  control  cards  for  consistency. 
Generates  some  arrays  and  files  based 
on  the  control  cards.  The  control  cards 
are  reformatted  and  output  on  file 
(ICARD)  for  final  processing  by  INPCRD 

Used  for  computing  interpolation  tables 
for  E.T.,  UTC  and  UTl  differences  (via 
call  to  subroutine  TIMARR) , and  computing 
and  sorting  time  intervals  for  which 
integrator  output  is  required.  Tables 
and  time  intervals  are  temporarily  stored 
in  file  32  (LUB) . (formerly  part  of 
subroutine  SPOLCD) 

Used  only  for  photogrammetric  data. 
Generates  ground  point  labels  for  output 
on  file  25  (IARCOT) , and  ground  point 
coordinates  output  on  direct  access  file 
30,  (130).  (entry  SP00L2)  clears  array 
for  ground  point  coordinates  (entry 
SP00L3).  (formerly  part  of  subroutine 
SPOLCD) 

Used  for  preliminary  processing  of  100 
series  cards  (entry  SP0100) . Used  for 
processing  meteorological  data  used  in 
refraction  formulae  (entry  SP0610)  and 
outputs  processed  data  as  well  as  TDRSS 
data  on  file  26  (IPASOT)  (entry  SP00L4) . 
(formerly  part  of  SPOLCD) 

Evaluate  integrated  power  series  output 
for  2nd  and  3rd  time  derivatives  of 
position 


1 


l 
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SUBROUTINE 


SUBROUTIN 

Name 

E 

CALLED 

BY 

PURPOSE  OR  FUNCTION 

SPOVEL 

SOFORT 

SOFSEC 

Evaluates  power  series  to  obtain 
satellite  state-vector 

STARAN 

PTSTAR 

Computes  predicted  stellar  camera 

orientation  angles  (photogrammetric  j 

measurement  types  1-6)  1 

! 

STARPA 

PTSTAR 

STARAN 

Computes  a (3x3)  matrix 
(See  note  1) 

STARPI 

STARAN 

Computes  a (3x3)  matrix 
(See  note  1) 

STASER 

TDRSS 

Develops  power  series  coefficients  for 
station  vector  in  inertial  space  at  a 
specified  time  T.  The  intertial  coordinate 
system  is  chosen  to  be  instantaneously 
coincident  with  the  "Earth  Fixed" 
coordinate  system  at  time  T 

STATEV 

ENGRAT 

Obtains  state-vector  (satellite  or 
planetary)  by  evaluating  previously 
developed  power  series  (*) 

SVAREQ 

SOFORT 

SOFSEC 

SOFRT 

Evaluates  power  series  to  obtain 
satellite  state  transition  matrix 

TDELAY 

TDRSS 

Computes  the  transmission  time  for  a 
radio  signal  sent  from  one  moving  point 
to  another 

TDRSS 

MEASXX 

Computes  predicted  TDRSS  measurements 
(Measurement  types  19-20)  j 

TERRA 

PTSTAR 

Computes  terrain  camera  orientation 
angles  from  stellar  camera  orientation 
angles 

TERRAS 

PTSTAR 

Computes  terrain  camera  orientation 
angles  such  that  the  camera  axes  point 
due  East  due  North  and  vertically  up 

TIMARR 

SPOOLl 

Rearranges  input  UT1  and  Ephemeris 

time  corrections  j 

UNIFD2 

EDIT 

Reads  data  tape  in  "unified"  Format  j 
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SUBROUTINE 


CALLED 

Name  BY  PURPOSE  OR  FUNCTION 


VISIBI 

SIMOUT 

VXPROD 

SELECT 

WRITER 

PTCMPA 

WTSPAX 

INP600 

INTCOD 

INPCRD 

XFORM 

TDRSS 

GEOCEG 

Computes  and  prints  time  of  satellite 
visibility  for  each  station 

Used  for  computing  the  vector  cross 
product  C = A * B. 

Utility  routine  for  writing  on  sequential 
file  (37) 

Utility  routine  for  direct  access 
file  41 


Uses  "Inertial  1950.0"  power  series  for 
satellite  state-vector  to  compute 
satellite  state-vector  power  series 
coefficients  in  same  coordinate  system 
as  used  by  STASER 


i 

| 


] 
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NOTES 


The  computed  matrices  are  of  the  form  of  matrices  (3.19.13) 
through  (3.9.16)  in  "A  Photogrammetric  and  Tracking  Network 
Analysis  Program,  Old  Dominion  Systems,  Inc.,  October  1973, 
Contract  DAAK  02-72-0-0434”. 

Routines  marked  "*"  are  used  in  the  integrator.  They  are 
called  once  per  integration  step.  Routines  marked  "**"  are 
used  in  the  integrator.  They  are  called  once  for  each  term 
(beyond  the  second)  in  the  power  series  expansion.  For  a 16 
term  power  series  expansion,  which  is  normal,  these  routines 
are  thus  called  14  times  per  integration  step. 


APPENDIX  G 


Photonap  Common  Blocks 

1 Listing  of  common  blocks  used  by  each  subroutine 

Routine  Common  Blocks 


ATMIN 


IATMOS 


ATMOUT 


ATMOS 


CATOOO 

DARITE 


DAYHMS 


TSPARM,  SOLFIL 


DBREAD 


TSPARM 


DBURN 

DEFALT 


DENS 


INTCMF,  INTCMO , CDEBUG , BURNS 

COMSOL,  ASPARM,  PSPARM,  GENCOM,  CONMET,  ICONST, 
EXTCM,  EARTH,  BURNS,  DRSGA,  SOLDRG 

IATMOS 


DGRITE 


FOTGNO 


DGRITS 


FOTGNO 


DRAGU 

DRAVAR 


DREDIT 


DUMDUM 


INTCMO,  DRSGB 

INTCMG , INTCMO,  INTCMI 

COMSOL,  ACINFO,  IONUMB,  GENCOM,  GPCOM,  CDEBUG, 
ICONST,  CWORK,  TSPARM 


EDIT 


EDITDT 


EIGN 

ENBVAR 

ENGRAT 


FOTGND,  COMSOL,  STINFO,  GENCOM,  EXTCM,  COVAR, 
CDEBUG,  CWORK,  IONUMB,  ICONST 

COMSOL,  TSEDIT,  ASPARM,  ACINFO,  PSPARM,  GPCOM, 
ATMOS,  DRSGA 


INTCMF,  INTCMO,  INTCMI,  GENCOM 

TIMING,  INTCMF,  INTCMO,  INTCMI,  POWER,  POTREC, 
CDEBUG,  EXTCM,  GENCOM,  IONUMB,  AJPL,  SEROUT 


I 


Routine 

Common  Blocks 

ENROOT 

EXPAND 

INTCMF,  INTCMO,  INTCMI,  AJPL,  POWER 

FIN ALP 

GENCOM,  TYLE,  TSPARM,  ASPARF,  PSPRMF,  CWORK, 

IONUMB , STINFO,  FOTGND,  EARTH,  COVAR 

FINDOG 

MISCOM,  INTCMF,  INTCMO,  INTCMI,  POTREC,  TSPARM, 
CDEBUG,  GENCOM 

FIREAD 

TSPARM 

FLREAD 

TSPARM 

PHOTO 

INROOT 

GENFIL 

COMSOL,  TSPARM,  TSEDIT,  ASPARM,  ACINFO,  PSPARM, 
STINFO,  GENCOM,  CDEBUG,  CWORK,  IONUMB,  ICONST, 
GPCOM,  TIMING,  INROOT,  AJPL,  EXTCM , POTREC, 

EARTH,  BURNS,  COVAR,  SOLDRG , DRSGA 

GENBUG 

COMSOL,  ASPARM,  PSPARM,  CWORK,  ICONST,  GPCOM, 
EXTCM,  BURNS,  SOLDRG 

GEOCEG 

CMEASR,  PARSOM,  EXTCM,  RSUMR,  EARTH,  PRTLB , PRTEMP 

GEOS 

STINFO,  IONUMB 

GETTAP 

CETBL2,  INTCMO,  CETBL9 , REC3 

GPRSM2 

IONUMB,  STINFO 

GROUND 

FOTO,  EARTH 

HOPRFT 

PRTLB,  CMEASR,  RSUMR,  XPNDR 

IEMCOL 

COMSOL 

IEMSET 

COMSOL 

IEMVAL 

COMSOL 

INPCRD 

TSPARM,  ASPARM,  CDEBUG,  IONUMB,  ICONST,  EXTCM, 
BURNS,  INPCMA 

INP100 

COMSOL,  GENCOM,  IONUMB,  ICONST,  GPCOM,  EXTCM, 
EARTH,  ATMOS,  INPCMA 

INP200 

COVAR,  ACINFO,  COMSOL,  ICONST,  POTREC,  EXTCM, 
BURNS,  INPCMA,  DRSGA 

-117-, 


Routine 


Common  Blocks 


INP300 

STINFO, 

FCONST,  INPCMA 

INP600 

TSPARM, 

CWORK, 

TSEDIT,  ASPARM,  ACINFO,  PSPARM,  COMSOL, 
POTREC,  GPCOM,  INROOT,  BURNS,  INPCMA.  DRSGA 

INP700 

COMSOL, 

GENCOM,  FCONST, 

INROOT,  INPCMA 

INPBUG 

TSPARM,  TSEDIT,  ASPARM, 
TYLE,  COMSOL,  ICONST 

ACINFO,  PSPARM, 

STINFO, 

INTCOD 

INROOT, 

TSPARM,  TSEDIT, 

CWORK,  GENCOM, 

GPCOM,  POTREC 

INTEG 

BURNS, 

CINTEG 

DRSGB,  CWORK,  INTCMF , INTCMO,  INTCMI,  EXTCM, 

INTERP 

TIMING 

INTGA 

POTREC, 
INTCMI , 

TIMING,  TSPARM, 
POWER 

EXTCM,  INTCMF, 

INTCMO, 

INTP1 

TIMING 

INTP2 

TIMING 

INTPRT 

INVERT 

SOLCOM, 

TSSOLV,  SOLFIL 

• 

INVPRE 

SOLFIL, 

TSPARM 

FOTGND,  SOLCOM, 

IONUMB,  CWORK, 

TSSOLV, 

INVRTB 

SOLCOM, 

TSSOLV,  SOLFIL, 

GENCOM,  CWORK, 

STINFO 

INVRTC 

SOLCOM, 

TSSOLV,  SOLFIL, 

GENCOM,  CWORK, 

STINFO 

INVSYM 

INVSYS 

JACHIA 

TSPARM 

JULDAY 

KEPLER 

KICKER 

SEROUT, 
EXTCM , 
INTCMO , 

CDEBUG , CINTEG,  CWORK,  CONMET,  GENCOM, 

IONUMB , PARTY,  POTREC,  POWER,  INTCMF, 

INTCMI,  CETBL1 

KMNCON 

EARTH, 

KMANI, 

AJPL,  CPSYST,  KMANI,  KMAN2,  KMAN3,  KMAN4, 
CINTEG,  EXTCM,  CWORK 

I 


liii  ,anri 


Routine 


KMNEVA 


KMNIDE 

KMNINI 


KMNINV 

KMNMP1 

KMNMP2 

KMNMP3 

KMNOUT 

KMNRAN 

KMNSIM 

KMNSM2 

MASCON 

MATINV 

MATZEV 

MATZEX 


MAT3T3 


MEASUR 


MEASXX 


NBDNEX 


Common  Blocks 


Routine 


Common  Blocks 


NBDPEX 

NEWJPL 

NUTION 

04DATE 

OCCULT 

OFDATE 

OLDMAN 

PAGE 

PFSOLV 

PFVARY 

PLATEC 


INTCMF,  INTCMO , 
INTCMF,  INTCMO, 
TIMING 
AJPL 

INTCMF,  INTCMO, 
AJPL 

GENCOM,  IONUMB , 
GENCOM 
EXTCM , POWER 
INTCMF,  INTCMO, 
FOTO 


INTCMI 

INTCMI,  CETBL1 , CETBL2 , CETBL4 


INTCMI 

OVRLAY 


INTCMI,  POWER 


POTSOL 

POTVAR 

PREINT 

PREPAR 

PRTBG1 

PRTBG2 


INTCMF,  INTCMO,  INTCMI,  AJPL 

CWORK,  CINTEG,  TSPARM , EARTH,  IONUMB,  EXTCM,  AJPL, 
DRSGB 

GENCOM,  EARTH,  TIMING,  CDEBUG , CWORK,  CINTEG,  BURNS, 
DRSGB,  PARTY,  POWER,  EXTCM,  IONUMB,  OVRLAY,  COVAR 

CWORK 

CWORK,  GENCOM,  STINFO,  IONR 


PRTBG3 


CINTEG,  RSUMR,  PRTLB 


PRTIAL 


CDEBUG,  CMEASR,  CONMET,  CWORK,  CINTEG,  BURNS, 
DRSGB,  PARTY,  SDP,  AJPL,  POTREC,  POWER,  GENCOM. 
EXTCM,  IONUMB,  STINFO,  TSPARM,  EARTH,  RSUMR, 

XPNDR,  IONR,  OVRLAY,  TIMING,  FOTGND , PRTLB,  PRTEMP 


PTCMPA 


PRTEMP,  STINFO,  CWORK,  GENCOM 


PTINIT 


STINFO,  PRTEMP,  CWORK,  CINTEG,  GENCOM,  EXTCM, 
TSPARM,  EARTH,  IONR,  RSUMR 


PTSTAR 


PRTEMP,  PRTLB,  CWORK,  CINTEG,  AJPL,  EXTCM,  IONUMB, 
EARTH,  TIMING,  FOTO,  FOTGND 
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Routine 


RAN601 

READE 

READER 

RECCOF 

RECPOT 

REDISK 

REDSKK 

REDS KM 

REFRCT 

REPRT2 

REREAD 

RESID 

ROOTDT 

ROTATE 

ROTINT 

ROTPAR 

ROTVFD 

ROT1 

RSUM 

SELECT 

SIGWT 

SIMOUT 

SKRIV 

SOFORT 


Common  Blocks 

CETBL2 , CETBL5 , CETBL1,  INTCMO,  CETBL4 , CETBL9 

INTCMF,  INTCMO,  INTCMI 

INTCMF,  INTCMO,  INTCMI,  CWORK,  POTREC,  GENCOM, 
AJPL 

TSPARM 

TSPARM,  TSSOLV 
TSSOLV,  SOLCOM,  SOLFIL 
EARTH,  GENCOM,  CMEASR,  CONMET 

TSPARM 

TSPARM,  CWORK,  STINFO,  GENCOM,  CDEBUG , IONUMB 

TIMING,  STINFO,  TYLE,  GENCOM,  CDEBUG,  MISCOM, 
CWORK,  IONUMB,  EXTCM,  EARTH,  ICONST,  POTREC, 
FCONST,  CONMET,  COVAR,  FOTGND,  TSPARM,  IONR 

OVRLAY,  STINFO,  EARTH 
EARTH,  PRTLB 
EARTH,  PRTLB 

RSUMR,  CMEASR,  PRTLB 
GPSYST,  KMANI,  KMAN2 

INROOT,  CWORK,  STINFO,  GENCOM,  IONUMB 

GENCOM,  CDEBUG,  PRTLB,  CINTEG,  EXTCM,  SDP,  POTREC 
BURNS,  IONUMB,  PARCOM,  POWER,  AJPL,  MISCOM 


Routine 


SOFRT 


SOFSEC 


SOLAWA 

SOLREC 

SOLVDT 

SOLVER 


SOLVFU 
SOL  VI 


SOREAD 

SORITE 

SPODT 

SPOLCD 


SPOOL 1 
SPOOL2 
SPOVEL 
SP0100 


SP02A3 

STARAN 

STARPA 


STARPI 


Common  Blocks 


GENCOM,  CDEBUG,  CINTEG,  EXTCM , IONUMB,  FINCOM, 
AJPL,  COVAR,  MISCOM 

GENCOM,  CDEBUG,  PRTLB , CINTEG,  EXTCM,  SDP,  POTREC, 
IONUMB,  PARCOM,  PARTY,  AJPL,  MISCOM 

CWORK 


SOLCOM,  TSSOLV,  SOLFIL 

SOLCOM,  FOTGND,  TSPARM,  TSSOLV,  IONUMB,  GENCOM, 
CWORK,  CDEBUG,  SOLFIL 

SOLCOM,  TSSOLV,  IONUMB,  GENCOM,  CWORK 

SOLCOM.  TSPARM,  TSSOLV,  GENCOM,  SOLFIL 

SOLCOM,  TSPARM,  TSSOLV,  GENCOM,  CWORK,  COVAR, 

SOLFIL 

GENCOM,  SOLFIL 
CWORK,  SOLFIL 
CWORK,  SOLFIL 
SP1COM,  SP2COM,  SP3COM 

CONMET,  TIMING,  INROOT,  FOTGND,  TSPARM,  STINFO, 
GENCOM,  IONUMB,  ICONST,  CWORK,  EARTH,  IONR, 
SP1COM,  SP2COM,  SP3COM 

TIMING,  IONUMB,  SP1COM 

FOTGND,  GENCOM,  IONUMB,  SP2COM 


CONMET,  POTREC,  INROOT,  TSPARM,  STINFO,  TYLE, 
GENCOM,  IONUMB,  CDEBUG,  CWORK,  SP3COM 


PRTLB,  PARCOM,  RSUMR,  EXTCM,  CMEASR,  EARTH,  PRTEMP 


Routine 


Common  Blocks 


STASER 

EARTH 

STATEV 

INTCMO , INTCMI 

SVAREQ 

TDELAY 

RSUMR 

TDRSS 

CMEASR,  PRTLB , PARCOM,  EXTCM , RSUMR,  XPNDR,  EARTH, 

GENCOM 

TERRA 

FOTO 

TERRAS 

FOTO,  PRTLB 

TIMARR 

TIMING 

UNIFD2 

STINFO,  IONUMB 

VISIBI 

GENCOM,  IONUMB,  STINFO 

VXPROD 

WRITER 

WTSPAX 

TSPARM,  TSEDIT,  IONUMB 

XFORM 

EXTCM 

IONOSF 

CMEASR,  PRTLB,  STINFO,  IONR,  IONTM 

ION 

CMEASR,  IONR 

BETA 

READU 

IONR 

DKSICO 

DKGK 

GK 

MAGFIN 

SICOJT 
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J 


G.2 


Listing  of  subroutines  utilizing  each  common  block 


Common 

Block 

Routines  Used  in 

ACINFO 

DREDIT,  EDITDT,  GENFIL,  INP200,  INP600,  INPBUG 

AJPL 

ENGRAT,  EXPAND,  GENFIL,  MEASUR,  MEASXX,  04DATE 
OFDATE,  POTVAR,  PREINT,  PRTIAL,  PTSTAR,  RECPOT, 
SOFORT,  SOFRT,  SOFSEC,  KMNCON 

ASPARF 

F INALP 

ASPARM 

DEFALT,  EDITDT,  GENFIL,  GENBUG,  INPCRD,  INP600, 
INPBUG 

ATMOS 

ATMOUT,  EDITDT,  INP100 

BURNS 

DBURN,  DEFALT,  GENFIL,  GENBUG,  INPCRD,  INP200, 
INP600,  INTEG,  PREPAR,  PRTIAL,  SOFORT 

CDEBUG 

DBURN,  DREDIT,  EDIT,  ENGRAT,  FINDOG,  GENFIL, 
INPCRD,  KICKER,  PREPAR,  PRTIAL,  RESID,  ROOTDT , 
SOFORT,  SOFRT,  SOFSEC,  SOLVER,  SP0100 

CETBL1 

KICKER,  NEWJPL,  READE 

CETBL2 

GETTAP,  NEWJPL,  READE 

CETBL4 

NEWJPL,  READE 

CETBL5 

READE 

CETBL9 

GETTAP,  READE 

CINTEG 

INTEG,  KICKER,  MEASUR,  MEASXX,  PREINT,  PREPAR, 
PRTBG3 , PRTIAL,  PTINIT,  PTSTAR,  SOFORT,  SOFRT, 
SOFSEC,  KMNCON 

CMEASR 

GEOCEG,  HOPRFT , MEASUR,  MEASXX,  PRTIAL,  REFRCT , 
RSUM,  SP02A3,  TDRSS,  IONOSF,  ION 

COMSOL 

DEFALT,  DREDIT,  EDIT,  EDITDT,  GENFIL,  GENBUG, 
IEMCOL , IEMSET,  IEMVAL,  INP100,  INP200,  INP600, 
INP700,  INPBUG 

CONMET 

DEFALT,  KICKER,  MEASUR,  PARTIAL,  REFRCT,  ROOTDT, 
SPOLCD,  SP0100 

COVAR 

EDIT,  FINALP,  GENFIL,  INP200,  PREPAR,  ROOTDT, 
SOFRT,  SOLV2 
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Common 

Block 

CWORK 


DRSGA 

DRSGB 

EARTH 

EXTCM 


FCONST 

FINCOM 

FOTGND 

FOTO 

GENCOM 


GPCOM 

GPSYST 

IATMOS 

ICONST 


Routines  Used  in ] 

SOLVFU,  SP0100,  DREDIT,  EDIT,  FINALP,  GENFIL,  J 

GENBUG,  INP600 , INTCOD,  INTEG,  INVPRE,  INVRTB , , 

INVRTC,  KICKER,  MASCON,  PREINT,  PREPAR,  PRTBG1 , 

PRTBG2 , PRTIAL,  PTCMPA,  PTINIT,  PTSTAR,  RECPOT, 

RES ID,  ROOTDT , SIMOUT,  SOLVER,  SOREAD,  SOLFIL, 

SPOLCD,  KMNCON,  KMNINI,  SOLAWA,  SOLV2 

DEFALT,  EDITDT,  INP200,  INP600,  GENFIL 

INTEG,  PREINT,  PREPAR,  PRTIAL,  DRAGU 

DEFALT,  FINALP,  GENFIL,  GEOCEG , GROUND,  INP100, 

MEASUR,  MEASXX,  PREINT,  PREPAR,  PRTIAL,  PTINIT, 

PTSTAR,  REFRCT,  ROOTDT,  ROTINT,  ROTPAR,  ROTVFD , 

SPOLCD,  SP02A3,  STASER,  TDRSS,  KMNCON 

DEFALT,  EDIT,  ENGRAT,  GENFIL,  GENBUG,  GEOCEG, 

INPCRD,  INP100,  INP200,  INTEG,  INTGA,  KICKER, 

MEASUR,  PFSOLV , PREINT,  PREPAR,  PRTIAL,  PTINIT, 

PTSTAR,  ROOTDT,  SOFORT,  SOFRT , SOFSEC,  SP02A3, 

TDRSS,  XFORM,  KMNCON,  KMNEVA  ' 

INP300,  INP700,  ROOTDT 


SOFRT 


DGRITE,  DGRITS,  EDIT,  FINALP,  INVPRE,  PRTIAL, 
PTSTAR,  ROOTDT,  SOLVER,  SPOLCD,  SPOOL2 

GROUND,  PLATEC,  PTSTAR,  STARAN,  TERRA,  TERRAS 

SOLV3 , SOLVFU,  SPOOL2 , SP0100,  VISIBI , KMNINI, 
KMNOUT , KMNSIM,  KMNSM2 , SOLV1,  SOLV2 , DEFALT, 
DREDIT,  EDIT,  ENBVAR,  ENGRAT,  FINALP,  FINDOG, 
GENFIL,  INP100,  INP700,  INTCOD,  INVRTB,  INVRTC, 
KICKER,  MEASUR,  MEASXX,  OLDMAN,  PAGE,  PREPAR, 
PRTBG2,  PRTIAL,  PTCMPA,  PTINIT,  RECPOT,  REFRCT, 
RES ID,  ROOTDT,  SIMOUT,  SOFORT,  SOFRT,  SOFSEC, 
SOLVER,  SPOLCD,  TDRSS 

DREDIT,  EDITDT,  GENFIL,  GENBUG,  INP100,  INP600, 
INTCOD 

KMNCON,  KMNINI,  KMNOUT,  KMNSIM,  SELECT 
ATMIN,  DENS 

DEFALT,  DREDIT,  EDIT,  GENFIL,  GENBUG,  INPCRD, 
INP100,  INP200 , INPBUG,  ROOTDT,  SPOLCD 
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Common 

Block 


Routines  Used  in 

INPCMA  INPCRD,  INP100,  INP200,  INP300,  INP600,  INP700 

INROOT  PHOTO,  GENFIL,  INP600,  INP700,  INTCOD,  SIMOUT, 

SPOLCD,  SP0100 

INTCMF  DBURN,  DRAVAR,  ENBVAR,  ENGRAT,  EXPAND,  FINDOG, 

INTEG,  INTGA,  KICKER,  MASCON,  MATZEX,  NBDNEX , 
NBDPEX,  NEW J PL,  OCCULT,  PFVARY , POTVAR,  RECCOF, 
RECPOT 

INTCMI  DRAVAR,  ENBVAR,  ENGRAT,  EXPAND,  FINDOG,  INTEG, 

INTGA,  KICKER,  MASCON,  MATZEV,  MATZEX,  NBDNEX, 
NBDPEX,  NEWJPL , OCCULT,  PFVARY,  POTVAR,  RECCOF, 
RECPOT,  STATE V 

INTCMO  DBURN,  DRAVAR,  ENBVAR,  ENGRAT,  EXPAND,  FINDOG, 

GETTAP,  INTEG,  INTGA,  KICKER,  MASCON,  MATZEV, 
MATZEX,  NBDNEX,  NBDPEX,  NEWJPL,  OCCULT,  PFVARY, 
POTVAR,  READE,  RECCOF,  RECPOT,  STATEV,  DRAGU 

IONR  PRTBG2 , PRTIAL,  PTINIT,  ROOTDT , SPOLCD,  IONOSF, 

ION,  READU 

IONTM  IONOSF 

IONUMB  KMNINI,  SOLVFU,  SPOOL1,  SPOOL2,  SPOIOO,  DREDIT, 

EDIT,  ENGRAT,  FINALP,  GENFIL,  GEOS,  GPRSM2 , 
INPCRD,  INP100,  INVPRE , KICKER,  OLDMAN,  PREINT, 
PREPAR,  PRTIAL,  PTSTAR,  RESID,  ROOTDT,  SIMOUT, 
SOFORT,  SOFRT,  SOFSEC,  SOLVER,  SPOLCD,  UNIFD2 , 
VISIBI , WTSPAX 


KMANI 

KMNCON,  KMNINI, 

KMNOUT, 

KMNSIM, 

KMNSM2 , SELECT 

KMAN1 

KMNCON,  KMNEVA, 

KMNINI, 

KMNINV 

KMAN2 

KMNCON,  KMNEVA, 

KMNINI, 

KMNOUT, 

KMNSM2 , SELECT 

KMAN3 

KMNCON 

KMAN4 

KMNCON,  KMNINI, 

KMNOUT, 

KMNSIM 

MISCOM 

FINDOG,  ROOTDT, 

SOFORT, 

SOFRT, 

SOFSEC 

OVRLAY 

OLDMAN,  PREPAR, 

PRTIAL, 

ROTINT 

PARCOM 

GEOCEG , SOFORT, 

SOFSEC, 

SP02A3 , 

TDRSS 

PARTY 

KICKER,  PREPAR, 

PRTIAL, 

SOFSEC 
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Block 

POTREC 

POWER 

PRTEMP 

PRTLB 

PSPARM 

PSPRMF 

REC3 

RSUMR 

SDP 

SEROUT 

SOLCOM 

SOLDRG 

SOLFIL 

SP1COM 

SP2COM 

SP3COM 

STINFO 

TIMING 


Routines  Used  in 


ENGRAT,  FINDOG,  GENFIL,  INP200,  INP600,  INTCOD, 
INTGA,  KICKER,  MASCON,  PRTIAL,  RECPOT,  ROOTDT , 
SOFORT,  SOFSEC,  SP0100 

ENGRAT,  EXPAND,  INTGA,  KICKER,  PFSOLV,  PFVARY , 
PREPAR,  PRTIAL,  SOFORT 

GEOCEG , MEASUR,  MEASXX,  PRTIAL,  PTCMPA,  PTINIT, 
FTSTAR , SP02A3 

GEOCEG,  HOPRFT,  MEASUR,  MEASXX,  PRTBG3 , PRTIAL, 
PTSTAR,  RGTPAR,  ROTVFD , RSUM,  SOFORT,  SOFSEC, 
SP02A3 , TDRSS , TERRAS,  IONOSF 

DEFALT,  EDITDT,  GENFIL,  GENBUG,  INP600,  INPBUG 

FINALP 

GETTAP 


GEOCEG,  HOPRFT,  MEASUR,  MEASXX,  PRTBG3 , PRTIAL, 
RSUM,  SP02A3,  TDELAY , TDRSS,  PTINIT 


PRTIAL, 

SOFORT, 

SOFSEC 

ENGRAT, 

KICKER 

INVERT, 

SOLVER, 

INVPRE , 
SOLV1 , 

INVRTB,  INVRTC,  REDSKM,  SOLVDT, 
SOLV2 , SOLVFU 

DEFALT, 

GENFIL, 

GENBUG 

DARITE , 
SOLVDT , 
SOLV3 

INVERT, 

SOLVER, 

INVPRE,  INVRTB,  INVRTC,  REDSKM, 
SOREAD,  SORITE , SOLV1,  SOLV2 , 

SPOLCD, 

SPODT, 

SPOOL1 

SPOLCD, 

SPODT, 

SPOOL2 

SPOLCD, 

SPODT, 

SP0100 

EDIT,  FINALP,  GENFIL,  GEOS,  GPRSM2 , INP300,  INPBUG, 
INVRTB , INVRTC , PRTBG2 , PRTIAL,  PTCMPA,  PTINIT, 
RESID,  ROOTDT,  ROTINT,  SIMOUT,  SPOLCD,  UNIFD2 , 
VISIBI , IONOSF,  KMNINI , SP0100 

ENGRAT,  GENFIL,  INTERI*,  INTGA,  INTP1,  INTP2 , 

NUTION , PREPAR,  PRTIAL,  PTSTAR,  ROOTDT,  SPOLCD, 
TIMARR,  SPOOL1 
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Common 

Block 


Routines  Used  in 


TSEDIT 

TSPARM 

TSSOLV 

TYLE 

XPNDR 


EDITDT,  GENFIL,  INP600 

DARITE , DBREAD,  DREDIT, 
FLREAD,  GENFIL,  INPCRD 
INTGA,  INVPRE , PREINT, 
REDSKK,  REREAD,  RES ID, 
WTSPAX,  JACHIA,  SOLV1, 

INVERT,  INVPRE,  INVRTB 
SOLVDT , SOLVER,  SOLVl, 

FINALP,  INPBUG,  ROOTDT 


INPBUG,  INTCOD,  WTSPAX 

FINALP,  FINDOC.  FIREAD, 
INP600,  INPBUG,  INTCOD, 
PRTIAL,  PTINIT,  REDISK, 
ROOTDT,  SOLVER,  SPOLCD, 
SOLV2 , SP0100 

INVRTC , REDSKK,  REDSKM, 
SOLV2 , SOLVFU 

SP0100 


HOPRFT , MEASXX,  PRTIAL,  TDRSS 
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